On mathematical problems of two-coefficient inverse problems of ultrasonic tomography

This article is concerned with developing effective methods for solving multi-coefficient inverse problems of wave tomography. Recently, there has been increased interest in tomographic imaging methods using ultrasonic, electromagnetic, and seismic radiation sources [1–5]. One interesting application of ultrasonic tomographic methods is early-stage breast cancer diagnosis [6–10]. Such applications require tomographic devices featuring high resolution and high accuracy. The most suitable mathematical model is the scalar wave model based on the hyperbolic differential equation, taking into account diffraction effects as well as absorption of ultrasonic waves. The inverse problem in this case is to reconstruct two unknown coefficients of the differential equation as functions of the coordinates.

Ultrasonic diagnostic methods are also widely used in nondestructive testing (NDT) [11–13]. A large number of studies on ultrasonic diagnostics in NDT and in medicine are focused on using reflected radiation. Ultrasonic tomographic methods being developed use the waves transmitted through the object in addition to reflected waves. Transmitted radiation provides the ability to perform quantitative characterization of inhomogeneities in the inspected objects with high accuracy. Such possibilities in real NDT experiments are presented in [14–16] assuming the scalar wave model. Some publications are concerned with solving inverse problems of tomographic ultrasonic diagnostics in a vector model, which supports transverse, surface, and other types of waves propagating in the medium in addition to the longitudinal wave [17–19].

This work is dedicated to the development of effective methods for solving two coefficient inverse problems in wave tomography under a model that takes into account both diffraction and absorption effects. From a mathematical point of view, the inverse problem is to minimize the residual functional between the wave field on the detectors, measured in the experiment, and the wave field computed using the given equation coefficients. The works [4, 20–26] provide representations for the gradient of the residual functional with mathematical rigor, which justifies minimizing the functional using iterative methods. Iterative methods of solving coefficient inverse problems for the wave equation based on minimizing the residual functional in application to seismic imaging are often referred to as the full waveform inversion (FWI) methods by the geophysical community [4]. The peculiarity of seismic imaging tasks is that typically seismic sources and receivers can be located only on the surface of the Earth. FWI approaches are also applicable to ultrasound imaging in NDT [11, 12].

The problem with wave tomography is that the functional is not convex and may have local minima. A large number of studies [16, 27–31] are dedicated to finding global minima of functionals. This phenomenon is referred to as cycle skipping in the FWI community. This problem is well known and various approaches have been developed to solve it. Data-domain methods are based on the modification of the misfit measurement [28]. Another approach uses multiple steps to solve the inverse problem. At the first stage, low-pass filtered signals are used, in subsequent stages the frequency is gradually increased. In various modifications, this approach is presented as the method of choice [30], multi-scale FWI in time or frequency domain [4, 32]. This paper uses this approach in the multistage method (MSM) as a gradient iterative method for solving wave tomography problems in the time domain. This paper provides estimates showing that at sufficiently low frequencies the nonlinear inverse problem approaches а linear problem and the residual functional approaches а convex one. This characteristic of inverse problems in wave tomography underlies the MSM method.

Another important problem in wave tomographic diagnostics is the uniqueness of the solution to inverse problems. The problem of uniqueness of solutions to coefficient inverse problems is explored in the works [33–38]. In this paper, a theorem is proven on the uniqueness of the solution to the inverse problem of ultrasonic tomography in ${{\boldsymbol{R}}^2}$ under a model that takes into account both diffraction effects and absorption. A reduction of the nonlinear inverse problem to a system of linear integral equations is proven for the limit approach as sounding frequency tends to zero.

Supercomputers are used to develop gradient iterative methods in this study. The effectiveness of the algorithms is illustrated by model calculations, which represent a case of tomographic diagnostics of soft tissues.

In this study, we consider waves described by the scalar wave equation with absorption taken into account. The scalar model facilitates computing the scalar wave field $u\left( {r,t} \right)\,$from given initial data using the equation

$$\begin{align}&\Delta u\left( {x,q,t} \right) - c\left( x \right){u_{tt}}\left( {x,q,t} \right) - a\left( x \right){u_t}\left( {x,q,t} \right) = - \delta \left( {x - q} \right)g\left( t \right),\end{align} \tag{ 1 }$$

$$\begin{align} \notag\\[-6pt] & u\left( {x,q,t = 0} \right) = {u_t}\left( {x,q,t = 0} \right) = 0.\end{align} \tag{ 2 }$$

Here ${c^{ - 0.5}}\left( x \right) = v\left( x \right)$ is the wave speed in the medium; $a\left( x \right)$ describes absorption in the medium; $x \in {{\boldsymbol{R}}^2}$; $\Delta$ is the Laplace operator on the variable $x$; $\delta$ is the Dirac function, which specifies the position of the point source at $q$. The sounding pulse generated by the source is described by the function.

The inverse problem of ultrasound tomography can be formulated as follows. Assume that the inhomogeneity of the medium is localized within the region $R$. Let us suppose that the positions of the wave sources, characterized by the parameter $q \in {{\boldsymbol{R}}^2}$, together form a set $Q$, and $Q$ does not intersect with $\bar R.$The measurements of the field $u\left( {x,q,t} \right)$ are available at points $d \in {{\boldsymbol{R}}^2}$, which collectively form a set $D$, which does not intersect with $\bar R$ (figure 1).

Zoom In Zoom Out Reset image size

We assume that $v\left( x \right)$, $a\left( x \right)$ are smooth functions in ${{\boldsymbol{R}}^2}$, $v\left( x \right)$ differs from the known constant ${v_0} = {\text{const}}$ only within the inhomogeneity of the medium ($0 \lt {v_1} \lt v\left( x \right) \lt {v_2}$), $a\left( r \right)$ differs from 0 also only within the inhomogeneity of the medium ($0 \unicode{x2A7D} a\left( x \right) \lt {a_2}$). For simplicity, we will set the function $g\left( t \right)$ that describes the sounding pulse to $g\left( t \right) = \delta \left( t \right)$. In the inverse problem, the objective is to determine the unknown functions $v\left( x \right)$ and $a\left( x \right)$ for $x \in R$, using the experimental data $U\left( {d,q,t} \right)$ obtained at points $D$ ($d \in D$) for various positions of the source $q \in Q$. The inverse problem of wave tomography in the given setting is a nonlinear coefficient inverse problem. The inverse problem of determining the coefficients of the wave equation was studied in works [39, 40].

We investigate the problem of uniqueness of the solution to the considered inverse problem in ${{\boldsymbol{R}}^2}$. To do this, we reduce the nonlinear inverse problem to an equivalent system of two uniquely solvable linear integral equations for the coefficients $c\left( x \right)$ and $a\left( x \right)$. This approach largely repeats the methods used in [33] for the wave equation containing not two, but only one unknown coefficient $c\left( x \right)$. The problem of the uniqueness of the solution to the inverse problem of determining two coefficients $c\left( x \right)$ and $a\left( x \right)$ in ${{\boldsymbol{R}}^3}$ was considered in [34], where the Laplace transform and the limit approach as $p \to 0$ for the Green's function of the operator ${c_0}{p^2}u\left( {x,p} \right) - \Delta u$ were applied to equation (1) with $p$ as the Laplace transform parameter. This technique is not suitable in the case of ${{\boldsymbol{R}}^2}$ space, since the Green's function for a two-dimensional space has a logarithmic singularity. The problem of the uniqueness of the solution to the inverse problem of determining a single coefficient $c\left( x \right)$ in ${{\boldsymbol{R}}^3}$ and the reduction of the nonlinear coefficient inverse problem to a linear integral equation of the first kind was considered in [35–38]. Let us present more rigorous mathematical formulations of the inverse problem.

Condition 1. Let $c\left( x \right)$ and $a\left( x \right)$ be twice continuously differentiable functions on ${{\boldsymbol{R}}^2}$. We also assume that the function $u\left( {x,t} \right)$ and its partial derivatives with respect to up to the second order have the Laplace transform defined for all $p \in {{\boldsymbol{C}}^ + }$, where ${{\boldsymbol{C}}^ + } = \{ p \in C:Re\,p > 0\}$, $\left( {{x_1},{x_2}} \right) = x \in {{\boldsymbol{R}}^2}$.

The Laplace transform is considered in the sense of generalized functions. As is known, this condition for the presence of the Laplace transform for the problem under consideration is not too limiting [33].

Condition 2. Let $R \subset {{\boldsymbol{R}}^2}$ be a bounded domain with a piecewise smooth boundary. D and Q are closed piecewise smooth curves that do not enclose the domain $R$ and do not intersect the closure $\bar R$.

We introduce the function $\xi \left( x \right) = \,c\left( x \right) - {c_0}\left( x \right)$, where ${c_0}^{ - 0.5}\left( x \right) = {v_0} = {\text{const}}$ is the wave velocity in a homogeneous medium surrounding the investigated inhomogeneity. Then, after the Laplace transform on equation (1), we obtain

$$\begin{align} \Delta \tilde u\left( {x,q,p} \right) - {p^2}\tilde u\left( {x,q,p} \right){c_0} = - \delta \left( {x - q} \right) + p\tilde u\left( {x,q,p} \right)a\left( x \right) + {p^2}\tilde u\left( {x,q,p} \right)\xi \left( x \right), \end{align} \tag{ 3 }$$

where $\tilde u\left( {x,q,p} \right)$ denotes the Laplace transform $\mathop {{{\mathop \int \nolimits}}}\nolimits_0^\infty \text{exp} \left( { - pt} \right)u\left( {x,q,t} \right){\text{d}}t$, $p \in {{\boldsymbol{C}}^ + }$. It is known that the Green's function for the operator on the left side of (3) has the following form in ${{\boldsymbol{R}}^2}$:

$$\begin{align}G\left( {x,q,p} \right) = - \frac{1}{{2\pi }}{K_0}\left( {\frac{{p\left| {x - q} \right|}}{{{v_0}}}} \right),\quad \quad \left| {\text{arg} p} \right| < \pi ,\end{align}$$

where ${K_0}\left( z \right)$ is the Macdonald function. Then equation (3) can be written as

$$\begin{align}\tilde u\left( {x,q,p} \right) + \frac{1}{{v_0^2}}G\left( {x,q,p} \right) = p\mathop \int \limits_R^{} G\left( {x,r,p} \right)\tilde u\left( {r,q,p} \right)\left( {a\left( r \right) + p\xi \left( r \right)} \right){\text{d}}r.\end{align} \tag{ 4 }$$

Dividing both sides of the equation (4) by $\text{ln} p$, where $p \in S_\varepsilon ^ + ,{ }S_\varepsilon ^ + = \left\{ p \in {\boldsymbol{C}}:\left| p \right| \le \varepsilon ,\right.$ $\left. Rep > 0 \right\},{ }$ $\varepsilon \lt 1$, we obtain the equation

$$\begin{align}\tilde v\left( {x,q,p} \right) + \frac{{G\left( {x,q,p} \right)}}{{v_0^2\text{ln} p}} = p\mathop \int \limits_{{S_\rho }}^{} G\left( {x,r,p} \right)\left( {a\left( r \right) + p\xi \left( r \right)} \right)\tilde v\left( {r,q,p} \right){\text{d}}r,\end{align} \tag{ 5 }$$

where $\tilde v\left( {x,q,p} \right) = {(\text{ln} p)^{ - 1}}\tilde u\left( {x,q,p} \right)$ and ${S_\rho }$ is a circle of radius $\rho$, such that $D,\,Q,\,R \subset {S_\rho }$.

In [33], the following representation was obtained for $p \to 0$:

$$\begin{align} G\left( {x,r,p} \right)& = - \frac{1}{{2\pi }}{K_0}\left( {\frac{{p\left| {x - r} \right|}}{{{v_0}}}} \right) \nonumber\\ & = \frac{1}{{2\pi }}\left( \text{ln} p + \gamma + \text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}} + \frac{{{{\left| {x - r} \right|}^2}}}{{4v_0^2}}{p^2}\text{ln} p + \frac{{{{\left| {x - r} \right|}^2}}}{{4v_0^2}}{p^2}\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}}\right. \nonumber\\ &\left. \quad + \frac{{\left( {\gamma - 1} \right){{\left| {x - r} \right|}^2}}}{{4v_0^2}}{p^2} \right) + O\left( {{{\left| p \right|}^3}} \right), \end{align} \tag{ 6 }$$

where $\gamma = 0,5772 \ldots$ is the Euler–Mascheroni constant.

Consider the operator ${B_p}:{L_2}\left( {{S_\rho }} \right) \to {L_2}\left( {{S_\rho }} \right)$

$$\begin{align}\left( {{B_p}w} \right)\left( x \right) = \mathop \int \limits_{{S_\rho }}^{} pG\left( {x,r,p} \right)\left( {a\left( r \right) + p\xi \left( r \right)} \right)w\left( r \right){\text{d}}r,\end{align} \tag{ 7 }$$

where $x \in {S_\rho }$, $p \in S_\varepsilon ^ +$.

Following the work [33], some properties of the operator ${B_p}$ can be clarified. Since $a\left( r \right) + p\xi \left( r \right)$ is continuous on ${S_\rho }$, and the kernel of the operator ${B_p}$ has the form (6), according to the properties of the potential-type operator from [41], it follows that the operator ${B_p}$ is completely continuous for all $p \in S_\varepsilon ^ + .$ From (6), it follows that $\left\| {{B_p}} \right\| = O\left( {\left| p \right|\left| {\text{ln} \left| p \right|} \right|} \right)$, so for a sufficiently small $\varepsilon$, for all $p \in S_\varepsilon ^ +$ the norm $\left\| {{B_p}} \right\| < 1$. Then for all $p \in S_\varepsilon ^ +$, denoting $g\left( {x,q,p} \right) = {( - v_0^2\text{ln} p)^{ - 1}}G\left( {x,q,p} \right)$, from (5) it follows

$$\begin{align}\tilde v\left( {x,q,p} \right) = {\left(I - {B_p}\right)^{ - 1}}g\left( {x,q,p} \right) = \left( {I + {B_p} + B_p^2 + \ldots } \right)g\left( {x,q,p} \right).\end{align} \tag{ 8 }$$

Due to the continuity of $G\left( {x,r,p} \right)$ by $p$, for simplicity we consider only real $p \in S_\varepsilon ^ +$. Then from (6) and (8), as well as using the representation

$$\begin{align}g\left( {x,q,p} \right) & = {g_0}\left( {x,q,p} \right) + O\left( {{{\left| p \right|}^3}{{\text{ln} }^{ - 1}}p} \right)\nonumber\\& = \frac{1}{{2\pi }}( 1 + \frac{\gamma }{{\text{ln} p}} + \frac{1}{{\text{ln} p}}\text{ln} \frac{{| {x - q} |}}{{2{v_0}}} + \frac{{{{| {x - q} |}^2}}}{{4v_0^2}}{p^2} + \frac{{{{| {x - q} |}^2}}}{{4v_0^2}}\frac{{{p^2}}}{{\text{ln} p}}\text{ln} \frac{{| {x - q} |}}{{2{v_0}}} \notag\\ & \quad + \frac{{( {\gamma - 1} ){{| {x - q} |}^2}}}{{4v_0^2}}\frac{{{p^2}}}{{\text{ln} p}} ) + O( {{{| p |}^3}{{\text{ln} }^{ - 1}}p} ).\end{align}$$

We can represent

$$\begin{align}\tilde v\left( {x,q,p} \right) & = {g_0}\left( {x,q,p} \right) + \frac{p}{{4{\pi ^2}}}{\left(1 + \frac{\gamma }{{\text{ln} p}}\right)^2}\text{ln} p\mathop \int \limits_R^{} \left( {a\left( r \right) + p\xi \left( r \right)} \right){\text{d}}r\notag\\ & \quad + \frac{p}{{4{\pi ^2}}}\left( {1 + \frac{\gamma }{{\text{ln} p}}} \right)\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}} + \text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)\left( {a\left( r \right) + p\xi \left( r \right)} \right){\text{d}}r\notag\\ & \quad + \frac{p}{{4{\pi ^2}\text{ln} p}}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}}\text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)\left( {a\left( r \right) + p\xi \left( r \right)} \right){\text{d}}r + O\left( {{{\left| p \right|}^3}{{\text{ln} }^{ - 1}}p} \right).\end{align} \tag{ 9 }$$

Assuming that $\varepsilon < {\text{exp}}\left( { - \gamma } \right)$, we can divide both sides by $\frac{p}{{4{\pi ^2}}}{(1 + \frac{\gamma }{{\text{ln} p}})^2}\text{ln} p$. As a result, similarly to [33], we have a power series expansion in terms of the small parameter ${(\text{ln} p + \gamma )^{ - 1}}$

$$\begin{align} h\left( {x,q,p} \right) & = \frac{{\tilde v\left( {x,q,p} \right) - {g_0}\left( {x,q,p} \right)}}{{\frac{p}{{4{\pi ^2}}}{{\left(1 + \frac{\gamma }{{\text{ln} p}}\right)}^2}\text{ln} p}}\notag\\ & = \mathop \int \limits_R^{} a\left( r \right) + p\xi \left( r \right){\text{d}}r + \frac{1}{{\text{ln} p + \gamma }}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}} + \text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)\left( {a\left( r \right) + p\xi \left( r \right)} \right){\text{d}}r\notag\\ & \quad + \frac{1}{{{{\left( {\text{ln} p + \gamma } \right)}^2}}}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}}\text{ln} \frac{{\left| {r - q} \right|}}{{2{c_0}}}} \right)\left( {a\left( r \right) + p\xi \left( r \right)} \right){\text{d}}r + O\left( {{{\left| p \right|}^2}} \right) = \mathop \int \limits_R^{} a\left( r \right){\text{d}}r\notag\\ & \quad + \frac{1}{{\text{ln} p + \gamma }}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}} + \text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)a\left( r \right){\text{d}}r \notag\\ & \quad + \frac{1}{{{{\left( {\text{ln} p + \gamma } \right)}^2}}}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}}\text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)a\left( r \right){\text{d}}r + O\left( {\left| p \right|} \right),\end{align} \tag{ 10 }$$

where the terms with $\xi \left( r \right)$ are included in $O\left( {\left| p \right|} \right).$ Evaluating the limits for $p \to 0$, we obtain a linear Fredholm integral equation of the first kind, which is called Lavrentiev's equation [35], for determining $a\left( r \right)$:

$$\begin{align}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}}\text{ln} \frac{{\left| {r - {\text{q}}} \right|}}{{2{v_0}}}} \right)a\left( r \right){\text{d}}r = {H_2}\left( {x,q} \right),\quad \quad q \in Q,\,\,x \in D,\end{align} \tag{ 11 }$$

where

$$\begin{align}{H_2}\left( {x,q} \right) = \mathop {\text{lim} }\limits_{p \to 0} {\left( {\text{ln} p + \gamma } \right)^2}\left( {h\left( {x,q,p} \right) - {H_0}\left( {x,q} \right) - \frac{{{H_1}\left( {x,q} \right)}}{{\text{ln} p + \gamma }}} \right);\end{align} \tag{ 12 }$$

$$\begin{align}{H_1}\left( {x,q} \right) = \mathop {\text{lim} }\limits_{p \to 0} \left( {\text{ln} p + \gamma } \right)\left( {h\left( {x,q,p} \right) - {H_0}\left( {x,q} \right)} \right);\quad \quad {H_0}\left( {x,q} \right) = \mathop {\text{lim} }\limits_{p \to 0} h\left( {x,q,p} \right);\end{align}$$

moreover, ${H_2}\left( {x,q} \right)$ can be computed from the data without knowing $a\left( x \right)$ and the limits exist. From equation (11) we find $a\left( x \right)$. Furthermore, the following equations are satisfied:

$$\begin{align}{H_1}\left( {x,q} \right) = \mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}} + \text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)a\left( r \right){\text{d}}r;\quad \quad {H_0}\left( {x,q} \right) = \mathop \int \limits_R^{} a\left( r \right){\text{d}}r.\end{align}$$

Next, $\xi \left( x \right)$ can be determined if $a\left( x \right)$ is known. Dividing equation (10) by $p$ and transferring the terms with $a\left( x \right)$ to the left side, we obtain:

$$\begin{align}m\left( {x,q,p} \right) & = \frac{1}{p}(h\left( {x,q,p} \right) - \mathop \int \limits_R^{} a\left( r \right){\text{d}}r - \frac{1}{{\text{ln} p + \gamma }}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}} + \text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)a\left( r \right){\text{d}}r\notag\\ & \quad - \frac{1}{{{{\left( {\text{ln} p + \gamma } \right)}^2}}}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}}\text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)a\left( r \right){\text{d}}r) = \mathop \int \limits_R^{} \xi \left( r \right){\text{d}}r\notag\\ & \quad + \frac{1}{{\text{ln} p + \gamma }}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}} + \text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)\xi \left( r \right){\text{d}}r \notag\\ & \quad + \frac{1}{{{{\left( {\text{ln} p + \gamma } \right)}^2}}}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}}\text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)\xi \left( r \right){\text{d}}r.\end{align} \tag{ 13 }$$

Similarly to (10), we have a power series expansion for $m\left( {x,q,p} \right)$ in terms of the small parameter ${(\text{ln} p + \gamma )^{ - 1}}$, where $q \in Q,\,\,x \in D,\,\,p \in \left( {0,\varepsilon } \right],\,\varepsilon \in \left( {0,{e^{ - \gamma }}} \right]$. Evaluating the limit for $p \to 0$, we obtain a linear Lavrentiev integral equation for determining $\xi \left( r \right)$:

$$\begin{align}\mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}}\text{ln} \frac{{\left| {r - {\text{q}}} \right|}}{{2{v_0}}}} \right)\xi \left( r \right){\text{d}}r = {M_2}\left( {x,q} \right),\,q \in Q,\,\,x \in D,\end{align} \tag{ 14 }$$

где

$$\begin{align} {M_2}\left( {x,q} \right) & = \mathop {\text{lim} }\limits_{p \to 0} {\left( {\text{ln} p + \gamma } \right)^2}\left( {m\left( {x,q,p} \right) - {M_0}\left( {x,q} \right) - \frac{{{M_1}\left( {x,q} \right)}}{{\text{ln} p + \gamma }}} \right)\notag\\ {M_1}\left( {x,q} \right) & = \mathop {\text{lim} }\limits_{p \to 0} \left( {\text{ln} p + \gamma } \right)\left( {m\left( {x,q,p} \right) - {M_0}\left( {x,q} \right)} \right);{M_0}\left( {x,q} \right) = \mathop {\text{lim} }\limits_{p \to 0} m\left( {x,q,p} \right),\end{align} \tag{ 15 }$$

where ${M_2}\left( {x,q} \right)$ depends on the previously obtained $a\left( x \right)$ but can be computed from the data without knowing $\xi \left( x \right),$ and the limits exist. Furthermore, the following equations are satisfied:

$$\begin{align}{M_1}\left( {x,q} \right) = \mathop \int \limits_R^{} \left( {\text{ln} \frac{{\left| {x - r} \right|}}{{2{v_0}}} + \text{ln} \frac{{\left| {r - q} \right|}}{{2{v_0}}}} \right)\xi \left( r \right){\text{d}}r;\quad {M_0}\left( {x,q} \right) = \mathop \int \limits_R^{} \xi \left( r \right){\text{d}}r.\end{align}$$

A significant number of works are concerned with proving the uniqueness of the solution of Lavrentiev's equation under various conditions imposed on the sets $Q,\,D$ in ${{\boldsymbol{R}}^2}$ and ${{\boldsymbol{R}}^3}$ [33, 37, 42]. These works provide proof of the theorem that equations of the type (11), (14) on ${{\boldsymbol{R}}^2}$ in the given statement have no more than one solution in the class C(R) of continuous functions in region R. Using this result, we formulate the following theorem.

Theorem 1. 

• (a)  
  Assume that conditions 1 and 2 for $a\left( x \right)$, $c\left( x \right)$, $u\left( {x,t} \right)$ and for sets $Q,\,D$, R from the problem statement are met. Then the inverse problem of finding the functions $a\left( x \right)$, $c\left( x \right)$ for equation (1) has a unique solution in the class C(R).
• (b)  
  The sought-for functions $a\left( x \right)$, $c\left( x \right)$ can be obtained via solving the linear integral equations (11), (12) and (14), (15) sequentially.

Theorem 1, besides proving the uniqueness of the solution to the inverse problem, also substantiates the reduction of the original nonlinear coefficient inverse problem to a system of two linear integral equations of the first kind with respect to functions $a\left( x \right)$ and $c\left( x \right)$. The numerical solution of linear integral equations is not a problem. However, the possibility of reducing the original nonlinear problem to a linear one depends on the presence of very low (close to 0) frequencies in the spectrum of the sounding signal. This is difficult to achieve in physical experiments. In addition, solving the obtained first-kind integral equations is a highly unstable problem due to the form of the kernel of the integral operator, among other things. The numerical solution of the inverse problem with integral operators of this type has a low spatial resolution that is unacceptable in tomographic imaging [43].

Iterative methods for minimizing the residual functional ${{\Phi }}\left( {c,a} \right)$ between the experimental data on the contour D and the wave field computed from equations (1) and (2) for given $c\left( x \right)$ and $a\left( x \right)$ have proven effective for solving inverse problems with respect to parameters $c\left( x \right)$ and $a\left( x \right)$ [44, 45].

$$\begin{align}{{\Phi }}\left( {c,a} \right) = \mathop \sum \limits_{j = 1}^M \frac{1}{2}\mathop \int \limits_0^T \mathop \int \limits_D^{} {\left( {u\left( {d,{q_{\boldsymbol{j}}},t;c,a} \right) - U\left( {d,{q_{\boldsymbol{j}}},t} \right)} \right)^2}{\text{d}}d{\text{ d}}t.\end{align} \tag{ 16 }$$

Here $U\left( {d,{q_{\boldsymbol{j}}},t} \right)$ are the experimental data on the contour $D$ over time (0,T), and $u\left( {x,{q_{\boldsymbol{j}}},t;c,a} \right)$ is the wave field obtained via solving the direct problems (1) and (2), which depends on the given coefficients $c\left( x \right)$ and $a\left( x \right)$. For multiple positions of the sounding radiation sources, the residual functional is a sum for $j = 1, \ldots ,M$ of the residual values obtained for each source position. For each fixed source $j$, the integral is summed over time (0,T) and over the contour $D$ containing all the detectors receiving the signals from the selected source. Mathematically, the inverse problem is posed as the problem of finding functions $\bar c\left( x \right)$ and $\bar a\left( x \right)$ that minimize the residual functional (16): $\bar c\left( x \right),\bar a\left( x \right):\mathop {\text{min} }\nolimits_{c\left( x \right),a\left( x \right)} {{\Phi }}\left( {c,a} \right) = {{\Phi }}\left( {\bar c,\bar a} \right)$. The functions $\bar c\left( x \right)$ and $\bar a\left( x \right)$ are taken as an approximate solution to the inverse problem.

The residual functional ${{\Phi }}\left( {c,a} \right)$ can be efficiently minimized using gradient methods. The gradient of the functional ${{\Phi }}\left( {c,a} \right)$ defined by expression (16) has the form [44, 45]:

$$\begin{align}{{\Phi }}_c^{\prime}\left( {c,a} \right) = \mathop \sum \limits_{j = 1}^M \mathop \int \limits_0^T \frac{1}{2}{u_t}\left( {x,{q_{\boldsymbol{j}}},t;c,a} \right){w_t}\left( {x,{q_{\boldsymbol{j}}},t;c,a} \right){\text{d}}t,\end{align}$$

$$\begin{align}{{\Phi }}_a^{\prime}\left( {c,a} \right) = \mathop \sum \limits_{j = 1}^M \mathop \int \limits_0^T \frac{1}{2}u\left( {x,{q_{\boldsymbol{j}}},t;c,a} \right){w_t}\left( {x,{q_{\boldsymbol{j}}},t;c,a} \right){\text{d}}t,\end{align} \tag{ 17 }$$

where $u\left( {x,{q_{\boldsymbol{j}}},t;c,a} \right)$ is the solution to the main problems (1) and (2), and $w\left( {x,{q_{\boldsymbol{j}}},t;c,a} \right)$ is the solution to the conjugate problems (18) and (19). Both solutions depend on coefficients $c\left( x \right)$ and $a\left( x \right)$.

$$\begin{align}& \Delta w\left( {x,{q_j},t} \right) - c\left( x \right){w_{tt}}\left( {x,{q_j},t} \right) + a\left( x \right){w_t}\left( {x,{q_j},t} \right) = \left( {u\left( {x,{q_j},t} \right) - U\left( {x,{q_j},t} \right)} \right){|_{x \in D}},\end{align} \tag{ 18 }$$

$$\begin{align}\notag\\[-6pt] & w\left( {x,{q_j},t = T} \right) = {w_t}\left( {x,{q_j},t = T} \right) = 0.\end{align} \tag{ 19 }$$

A typical situation for nonlinear problems is a non-convex residual functional ${{\Phi }}\left( {c,a} \right)$ that may have local minima. As a result, using gradient methods of functional minimization with arbitrary initial approximation guarantees the convergence only to a local minimum, but not to the global minimum.

Let us consider another approach to solving the problem of nonlinearity, not related to the reduction of the original nonlinear inverse problem to a system of linear integral equations of the first kind (11), (12) and (14), (15). Unlike the reduction to a linear problem, we will not evaluate the limits for $p \to 0$, but instead we will use sufficiently low (but not equal to 0) frequencies at the first stage. This constitutes the so-called MSM.

Let slightly stronger conditions than condition 1 be satisfied. We assume that the function $u\left( {x,t} \right)$ and its partial derivatives with respect to ${x_1},{x_2},\,t$ up to the second order have a Laplace transform defined for all $p \in {\boldsymbol{C}}_0^ +$, where ${\boldsymbol{C}}_0^ + = \left\{ {p \in C:Re\,p \ge 0,\,p \ne 0} \right\}$, $\left( {{x_1},{x_2}} \right) = x \in {{\boldsymbol{R}}^2}$. That is, it is additionally assumed that there is a Laplace transform for imaginary$\,p = iw$. The Laplace transform here is defined over generalized functions. Applicability of this condition was studied in detail in [33]. It was noted that this condition is not too restrictive. In particular, it holds for the Green's function $E\left( {x,t,q} \right)$ of the wave equation in homogeneous space

$$\begin{align}E\left( {x,t,q} \right) = \frac{{\theta \left( {t - \left| {x - q} \right|/{v_0}} \right)}}{{2\pi v_0^2\sqrt {{t^2} - {{\left(\left| {x - q} \right|/{v_0}\right)}^2}} }}\ .\end{align}$$

Using the approach described above in deriving equation (3), the Helmholtz equation with absorption can be represented as follows:

$$\begin{align}\Delta \tilde u\left( {x,q,w} \right) + \frac{{{w^2}}}{{v_0^2}}\tilde u\left( {x,q,w} \right) = - \delta \left( {x - q} \right) + iw\tilde u\left( {x,q,w} \right)a\left( x \right) - {w^2}\tilde u\left( {x,q,w} \right)\xi \left( x \right),\end{align} \tag{ 20 }$$

where $w$ is the sounding frequency, $\xi \left( x \right) = \,{v^{ - 2}}\left( x \right) - v_0^{ - 2}$. Using the Green's function method for the operator on the left side of equation (20), we can write:

$$\begin{align}\tilde u\left( {x,q,w} \right) + G\left( {x,q,w} \right) = iw\mathop \int \limits_R G\left( {x,r,w} \right)\tilde u\left( {r,q,w} \right)a\left( r \right){\text{d}}r - {w^2}\mathop \int \limits_R G\left( {x,r,w} \right)\tilde u\left( {r,q,w} \right)\xi \left( r \right){\text{d}}r,\end{align} \tag{ 21 }$$

where $G\left( {x,r,w} \right) = - \frac{i}{4}H_0^{\left( 2 \right)}\left( {k\left| {x - r} \right|} \right)$ is the Hankel function of the second kind, $k = \frac{w}{{{v_0}}}$. As is known, ${K_0}\left( {iw} \right) = - \frac{{\pi i}}{2}H_0^{\left( 2 \right)}\left( w \right)$ and ${K_0}\left( { - iw} \right) = \frac{{\pi i}}{2}H_0^{\left( 1 \right)}\left( w \right)$, where ${K_0}\left( z \right)$ is the Macdonald function from (6). Consequently, the estimates obtained earlier for the Macdonald function are valid for the Hankel functions.

Equation (21) is a nonlinear integral equation with respect to $a\left( r \right)$ and $\xi \left( r \right)$. Denoting $\zeta \left( r \right) = w\xi \left( r \right) - ia\left( r \right)$, the equation (21) can be represented in another form, which highlights the first term of the Born series.

$$\begin{align} & \tilde u\left( {x,q,w} \right) + G\left( {x,q,w} \right)\end{align}$$

$$\begin{align} & \quad = w\mathop \int \limits_R G\left( {x,r,w} \right)\zeta \left( r \right)G\left( {r,q,w} \right)dr + {w^2}\mathop \int \limits_R G\left( {x,r,w} \right)\zeta \left( r \right) \notag\\ & \qquad \times \mathop \int \limits_R G\left( {r,\gamma ,w} \right)\zeta \left( \gamma \right)\tilde u\left( {\gamma ,q,w} \right)d\gamma dr\,.\end{align} \tag{ 22 }$$

The first term on the right-hand side of (22) is the Born approximation of the scattering problem and takes the form of a linear integral operator with respect to $\zeta \left( r \right)$. The second term on the right-hand side of (22) has a multiplier ${w^2}$ and two multipliers $\zeta \left( r \right)\,$ under the integral. For sufficiently small values of $\zeta \left( r \right)\,$ and for $w \to 0$, the second term on the right-hand side of equation (22) becomes a small quantity compared to the first term, and equation (22) becomes close to a linear equation. Therefore, the functional (16) becomes close to a convex one.

Let us formulate this more rigorously by writing equation (22) as follows:

$$\begin{align}\begin{gathered} \tilde u\left( {x,q,w} \right) + G\left( {x,q,w} \right) \hfill \\ \quad = w\mathop \int \limits_R G\left( {x,r,w} \right)\zeta \left( r \right)\left( {G\left( {r,q,w} \right) + w\mathop \int \limits_R G\left( {r,\gamma ,w} \right)\tilde u\left( {\gamma ,q,w} \right)\zeta \left( \gamma \right){\text{d}}\gamma } \right){\text{d}}r \hfill \\ \end{gathered} .\end{align} \tag{ 23 }$$

Let us show that the following limit is zero:

$$\begin{align}\mathop {\text{lim} }\limits_{w \to 0} \frac{{\left\| {w\mathop \int \nolimits_R G\left( {r,\gamma ,w} \right)\tilde u\left( {\gamma ,q,w} \right)\zeta \left( \gamma \right)d\gamma } \right\|}}{{\left\| {G\left( {r,q,w} \right)} \right\|}} = 0\,,\end{align}$$

that is, for sufficiently small w, the inequality which we will call the condition of closeness to the linear problem holds:

$$\begin{align}\left\| {G\left( {r,q,w} \right)} \right\| \gg \,\left\| {w\mathop \int \limits_R G\left( {r,\gamma ,w} \right)\tilde u\left( {\gamma ,q,w} \right)\zeta \left( \gamma \right){\text{d}}\gamma } \right\|,\end{align} \tag{ 24 }$$

where ||⋅|| is the norm in ${L_2}\left( R \right)$. For the Green function $G\left( {r,q,w} \right)$, from (6) it follows that for all sufficiently small $w$ $\left\| {G\left( {r,q,w} \right)} \right\| > {m_3}$ for some ${m_3} > 0$. Due to the boundedness and continuity of functions $\xi \left( \gamma \right)$ and $a\left( \gamma \right)$, the inequality $\left| {\zeta \left( \gamma \right)} \right| < {m_1}$ holds for some ${m_1} > 0$. From the energy inequality for the Helmholtz equation, $\left\| {\tilde u\left( {\gamma ,q,w} \right)} \right\| < {m_2}$ for some ${m_2} > 0$. Consider the operator ${B_w}:{L_2}\left( R \right) \to {L_2}\left( {R\,} \right)$

$$\begin{align}\left( {{B_w}\chi } \right)\left( x \right) = - {w^{0.25}}\int\limits_R G\left( {x,r,w} \right)\chi \left( r \right){\text{d}}r\ ,\end{align}$$

where $r \in R$, $x \in R$. According to the properties of a potential-type operator from [41], the operator ${B_w}$ is completely continuous. From (6), it follows that $\left\| {{B_w}} \right\| = O( {{{\left| w \right|}^{0.25}}\left| {\text{ln} \left| w \right|} \right|} )$, so that $\left\| {{B_w}} \right\| < 1\,$for all sufficiently small $w$. Thus, for sufficiently small $w$ the following inequality holds:

$$\begin{align} & \left\| {w\mathop \int \limits_R G\left( {r,\gamma ,w} \right)\tilde u\left( {\gamma ,q,w} \right)\zeta \left( \gamma \right){\text{d}}\gamma } \right\| \unicode{x2A7D} {w^{0.75}}\left\| {{B_w}} \right\|{\left( {\mathop \int \limits_R {{\left( {\tilde u\left( {\gamma ,q,w} \right)\zeta \left( \gamma \right)} \right)}^2}{\text{d}}\gamma } \right)^{0.5}}\end{align}$$

$$\begin{align} & \quad \unicode{x2A7D} {w^{0.75}}{m_1}{\left( {\mathop \int \limits_R^{} {{\tilde u}^2}d\gamma } \right)^{0.5}} = {w^{0.75}}{m_1}\left\| {\tilde u\left( {\gamma ,q,w} \right)} \right\| \unicode{x2A7D} {w^{0.75}}{m_2}{m_1} \notag\\ & \qquad \ll {m_3} \lt \left\| {G\left( {r,q,w} \right)} \right\|.\end{align} \tag{ 25 }$$

The idea that for small heterogeneities $\zeta \left( r \right)$ and sufficiently small $w$ the nonlinear inverse problem approaches a linear problem and the residual functional becomes convex is at the core of the MSM. The MSM algorithm consists of the following steps. The initial approximation is chosen as constants $v\left( r \right) = {v_0}$ for the speed of sound and $a\left( r \right) = {a_0} = 0$ for the absorption factor, where ${v_0}$ is a known speed of sound in the homogeneous medium surrounding the imaged object. Experimental data at two (or more) different central frequencies ${w_1}$ and ${w_2}$ are used, ${w_1} < {w_2}$. First, the inverse problem is solved using iterative methods using the data for the lower frequency ${w_1}$, and a constant initial approximation. At low frequencies, the residual functional approaches convex and the iterative minimization process does not encounter local minima. Then, at the second stage, the inverse problem is solved using the higher frequency ${w_2}$ to improve the resolution. The speed of sound and absorption factor obtained via solving the inverse problem at the first stage using the lower frequency ${w_1}\,$are used as the initial approximation for the second stage. The process continues further if more than two central frequencies are chosen.

Note that in this article, experimental data with a central frequency $w$ is understood as data obtained by sounding with a broadband pulse with a central frequency $w$ and the problem is solved in the time domain. It is also possible to perform sounding with harmonic signals with different frequencies and solve the problem in the frequency domain for the Helmholtz equation.

In the provided estimates of the approximation of the nonlinear inverse problem with a linear one at sufficiently low frequencies, nothing is said about the specific values of the parameters of scattering objects and the frequencies in a specific problem. Let us try to estimate the scattering wave function depending on the parameters of the scattering object and the frequency, in order to determine at what values of $w$ and $\left( {ia\left( r \right) - w\xi \left( r \right)} \right)$ inequality (24) holds.

The integral in inequality (24) represents the scattered wave $\tilde u\left( r \right)$ in the area of heterogeneity $R \subset {{\boldsymbol{R}}^2}$. The expression for the scattered wave can be obtained by following the results of [46] for perturbation theory in quantum mechanics in the WKB approximation. Let us consider the equation for the scattering wave function $\tilde u\left( r \right)$ for a small scatterer $\zeta \left( r \right)$, $\,r = \left( {x,y} \right) \in {{\boldsymbol{R}}^2}$, which is similar to equation (20) without taking into account the point source of radiation

$$\begin{align}\Delta \tilde u\left( r \right) + {k^2}\tilde u\left( r \right) = - \tilde u\left( r \right)w\zeta \left( r \right),\end{align} \tag{ 26 }$$

where $k = \frac{w}{{{v_0}}}$. Let the incident plane wave propagate along the х-axis. In the first approximation, the wave function depends on the coordinates in the same way as the incident plane wave propagating along the x-axis. Therefore, we seek a solution in the form $\tilde u\left( r \right) = F\left( r \right)\,{\text{exp}}\left( { - ikx} \right)$, where $F\left( r \right)$ changes slowly compared to ${\text{exp}}\left( { - ikx} \right)$. Substituting $\tilde u\left( r \right)$ into equation (26), we obtain:

$$\begin{align} & {F_{xx}}\text{exp} \left( { - ikx} \right) - 2ik{F_x}\text{exp} \left( { - ikx} \right) - {k^2}F\text{exp} \left( {ikx} \right) + {F_{yy}}\text{exp} \left( { - ikx} \right) + {k^2}F\text{exp} \left( { - ikx} \right)\notag\\ & \quad = - w\zeta \left( {\boldsymbol{x}} \right)F\left( r \right)\,{\text{exp}}\left( { - ikx} \right)\ .\end{align}$$

We simplify this expression by dividing by $\text{exp} \left( { - ikx} \right).$ For large $k$ and slowly varying $\,F\left( r \right)$ we can neglect the second derivatives of $F\left( r \right)$, then

$$\begin{align} - 2ik{F_x}\left( r \right) = - w\zeta \left( r \right)F\left( r \right);\quad \quad F\left( {x,y} \right) = {\text{exp}} - \frac{{iw}}{{2k}}\int\limits_{ - \infty }^x \zeta \left( {\chi ,y} \right){\text{d}}\chi .\end{align}$$

So,

$$\begin{align}\tilde u\left( r \right) = F\left( r \right)\text{exp} \left( { - ikx} \right) = \text{exp} \left( { - i(kx + \frac{w}{{2k}}\mathop \int \limits_{ - \infty }^x \zeta \left( {\chi ,y} \right)d\chi } \right).\end{align} \tag{ 27 }$$

In physical media heterogeneous in speed and absorption, according to formula (21), the scattering function has the form $\zeta \left( r \right) = - ia\left( r \right) + w\xi \left( r \right)$. Substituting $\zeta \left( r \right)$ into (27), we get

$$\begin{align}\tilde u\left( r \right) & = \text{exp} \left( { - i(kx + \frac{1}{{2k}}\mathop \int \limits_{ - \infty }^x - iwa\left( {\chi ,y} \right) + {w^2}\xi \left( {\chi ,y} \right){\text{d}}\chi } \right)\end{align}$$

$$\begin{align} & \approx \text{exp} \left( { - i\left( {kx + \frac{{{v_0}}}{2}w\tilde \xi \tilde R} \right) - \frac{{{v_0}}}{2}\tilde a\tilde R} \right),\end{align} \tag{ 28 }$$

where $\tilde R$ is the size of the heterogeneity $R$ from the beginning of the heterogeneity to point $r$ along the line parallel to the x-axis and passing through $r$; $\tilde \xi$ and $\tilde a$ are average values of parameters over the heterogeneity. Thus, according to (28), heterogeneity in speed leads to an additional phase accumulation of $\frac{{{v_0}}}{2}w\tilde \xi \tilde R$, and heterogeneity in absorption leads to a decrease in amplitude with a factor of $\text{exp} \left( { - \frac{{{v_0}}}{2}\tilde a\tilde R} \right)$.

To meet the condition of proximity to the linear problem (24) for the incident plane wave $\text{exp} \left( { - ikx} \right)$ in the region of heterogeneity $R$, the following inequality must hold:

$$\begin{align}\left\| {\tilde u\left( r \right) - \text{exp} \left( { - ikx} \right)} \right\| \ll \left\| {\text{exp} \left( { - ikx} \right)} \right\|.\end{align} \tag{ 29 }$$

Thus,

$$\begin{align} & \tilde u\left( r \right) - \text{exp} \left( { - ikx} \right) = \text{exp} \left( { - i\left( {kx + \frac{{{v_0}}}{2}w\tilde \xi \tilde R} \right) - \frac{{{v_0}}}{2}\tilde a\tilde R} \right) - \text{exp} \left( { - ikx} \right)\end{align}$$

$$\begin{align} & \quad \approx \text{exp} \left( { - ikx} \right)\left( {1 - i\frac{{{v_0}}}{2}w\tilde \xi \tilde R - \frac{{{v_0}}}{2}\tilde a\tilde R - 1} \right)\ .\end{align}$$

Thus, when the condition is met

$$\begin{align}\left| {i\frac{{{v_0}}}{2}w\tilde \xi \tilde R + \frac{{{v_0}}}{2}\tilde a\tilde R} \right| \ll 1\end{align} \tag{ 30 }$$

condition (29) is also fulfilled. Note that the condition $\left| {\frac{{{v_0}}}{2}w\tilde \xi \tilde R} \right| \ll 1$ coincides with the well-known conditions for the phase accumulation in wave problems [16, 27–31].

Conditions similar to inequalities (24) and (30) in quantum mechanics, as well as in [47] are called the conditions of applicability of the Born approximation. That is, when these conditions are met, it is proposed to solve the inverse problem in the Born approximation. The fundamental difference of the approach considered in this study is that the Born approximate formulation is not used to solve the inverse problem. At all stages, the problem is solved in a full formulation via wave propagation simulations. But at the first stage, a sufficiently low sounding frequency is chosen, at which the condition of proximity to the linear problem (24) is met. At subsequent stages, condition (24) can be relaxed by applying higher frequencies, since an initial approximation for these stages becomes close to the exact solution.

The capabilities of the MSM method are demonstrated on model problems. The finite-difference time-domain method on a uniform rectangular grid in ${{\boldsymbol{R}}^2}$ was used to solve equations (1) and (2). The finite difference scheme approximates equation (1) with second order accuracy.

The values of the parameters of model problems were chosen typical for the problem of breast cancer diagnosis. At the first stage of the MSM method, a sounding pulse with an average wavelength of 15 mm and an initial approximation of $c\left( x \right) = {c_0} = {\text{const}}$ and $a\left( x \right) = 0$ were used. The wavelength at the first stage should be large enough to ensure the convergence of the gradient descent method started from the chosen initial approximation. At the next stage, a pulse with an average wavelength of 5 mm was used, and the initial approximation was the approximate solution obtained at the first stage.

Figure 2 shows the scheme of the numerical tomographic experiment. 24 sources of sounding waves were evenly spaced at the boundary of the square computational region $L$. Ultrasonic wave detectors were also located on the boundary of region $L$ with a step less than 1 mm. An object $R$ with a speed of sound and absorption differing from the background is located in region $L$ and is surrounded by sources and detectors. The size of the object $R$ is about 10 cm, the side of the square region $L$ is 15 cm. The speed of sound inside the object differs from the speed in the surrounding medium ${v_0} = 1.5$ km s⁻¹ by no more than 10%.

Zoom In Zoom Out Reset image size

Remark 1. In the numerical experiment, sources and receivers are located at a set of points on the boundary of the region $L$ that encloses $R$. This arrangement of sources and receivers on curves enclosing $R$ is common in real experiments. However, the location of the boundary of domain $L$ differs from condition 2 of the uniqueness theorem 1, where the sources and receivers are located on curves $D$ and $Q$ that do not enclose the domain $R$. The proof of theorem 1 uses the uniqueness theorem for the solution of Lavrentiev's equation under the assumption that curves $D$ and $Q$ do not enclose $R$. The authors are not aware of any proof of the uniqueness of the solution of Lavrentiev's equation in ${{\boldsymbol{R}}^2}$ under the assumption that curves $D$ and $Q$ enclose $R$. The technique of proving the uniqueness of the solution of Lavrentiev's equation under condition 2 is not suitable in the case when $D$ and $Q$ enclose $R$, since the following is important in the proof: let us denote a logarithmic potential with support represented by a harmonic function in $R$ as $P\left( x \right)$. Let $P\left( x \right) = 0$ for $x \in D$ or $x \in Q$, then the proof requires that $P\left( x \right) = 0$ for $x \notin R$ [48]. This is true when $D$ and $Q$ do not enclose the domain $R$. When $D$ and $Q$ enclose the domain, this may not be true, due to the exterior boundary value problem in ${{\boldsymbol{R}}^2}$. However, from numerics it is clear that, most likely, the uniqueness theorem is also satisfied if $D$ and $Q$ enclose the domain $R$. We are convinced that the reconstruction results are much better when the sources and receivers are located around $R$. This can also be seen from comparisons of figures 4(c) and 6. At the same time, it is obvious that the regions $D$ and $Q$ can be curved so that they actually surround $R$, look like narrow strips and do not enclose the domain $R$.

Figure 3(a) shows the sounding pulses used for calculations with average wavelengths of 5 mm and 15 mm. Figure 3(b) shows the frequency spectra of these pulses. In physical tomographic experiments, sounding can be conducted either with different emitters for wavelengths of 5 mm and 15 mm or with a wideband pulses with an average wavelength of 5 mm, and then the higher frequencies can be filtered out to obtain the signal with a wavelength of 15 mm.

Zoom In Zoom Out Reset image size

To conduct model experiments, a computational grid of 702 × 702 points with a step of 0.21 mm was chosen. Waves were measured by receivers at each grid point along the perimeter of square $L$. The sources have a wide radiation pattern, namely, the sources are located on the boundary of square $L$ and emit a pulse in the form of a semicircle into the area $L$. The sources take turns emitting a sounding pulse, which is registered by all of the receivers. The registration time T from formula (16) was chosen to ensure registration of the pulse from any source by all of the receivers. At the boundary of region $L$, the condition of transparency of the boundaries was applied. The results of solving the direct problem for each source, namely the data registered by each receiver, are stored and used as experimental data.

Explicit expression (17) for the gradient of the residual functional was used to solve the inverse problem via an iterative gradient-descent method. A priori information that near the sides of the square $L$ the velocity and absorption factor are equal to $c\left( x \right) = {c_0}$ and $a\left( x \right) = 0$, respectively, was used, so after each iteration the values near the boundary were smoothly clamped to these values. The iterative process ends when the output image stops changing. Typically, about 100–200 iterations were used at each stage of the method. Reconstructing only the velocity requires 2–3 times fewer iterations; the absorption factor is reconstructed much more slowly. Calculations were carried out on 100 computing cores of Intel Haswell-EP, 2.6 GHz CPUs on 'Lomonosov-2' supercomputer at Lomonosov Moscow State University. The calculation time for one inverse problem was about 1 h.

Figure 4 shows phantoms (figures 4(a) and (d)) and the reconstructed images of sound speed and absorption factor obtained using the MSM method. Two stages of the method were used for image reconstruction. At the first stage, the spectrum of the sounding signal was chosen with a central frequency of 100 kHz and an average wavelength of 15 mm, which ensures the convergence of the iterative gradient descent process. The reconstruction results for the first stage are shown in figures 4(b) and (e). Using these images as an initial approximation, the second stage of the iterative process was performed with an average sounding frequency of 300 kHz, corresponding to a wavelength of 5 mm. The reconstruction results for the second stage are shown in figures 4(c) and (f). These reconstruction results show that the first stage ensured the convergence of the iterative process to some approximate solution, however, the spatial resolution at the first stage is significantly lower than the spatial resolution obtained at the second stage. In figure 4(c), four point-like objects 3.5 mm in diameter are clearly visible. As can be seen from the comparison of figures 4(c) and (f), the reconstruction accuracy of the absorption coefficient is inferior to that of the speed of sound. This can be explained by the fact that the absolute value of the gradient with respect to the speed of sound in formula (5) is greater than that of the gradient with respect to the absorption factor.

Zoom In Zoom Out Reset image size

Figure 5(a) shows one-dimensional graphs of the residual functional (4) along the line passing through the initial approximation at $\alpha = 0$ and the exact solution of the problem at $\alpha = 1\,$ for the average wavelength $\lambda = 15\,$ mm (dashed line) and $\,\,\lambda = 5\,$ mm (solid line) depending on the parameter $\alpha$. The simulation parameters for each $\alpha$ are calculated as $c\left( {x;\alpha } \right) = {c_0}\left( {1 - \alpha } \right) + \bar c\left( x \right)\alpha$, $a\left( {x;\alpha } \right) = {a_0}\left( {1 - \alpha } \right) + \bar a\left( x \right)\alpha$, where ${c_0},\,{a_0}$ are parameters in water, and $\bar c\left( x \right)$ and $\bar a\left( x \right)$ constitute the exact solution of the problem. Since the calculations were carried out with exact simulated data without measurement errors, the residual functional is zero at the exact solution of the inverse problem. As can be seen from the graph in figure 5(a), the residual functional plot for the wavelength of 5 mm has a bend, which may be associated with a local minimum. Figure 5(a) indirectly shows that the longer the wavelength, the wider the area of convergence of iterative gradient algorithms to the approximate solution.

Zoom In Zoom Out Reset image size

Figures 5(b) and (c) show the reconstruction results of the speed of sound and absorption coefficient obtained via the iterative gradient method using a wideband signal with an average wavelength of 5 mm and an initial approximation ${c_0},\,{a_0}$. Given a short wavelength of the sounding pulse, the iterative process started from ${c_0},\,{a_0}$ stops at the local minimum of the residual functional, resulting in an incorrect solution.

The reconstructed images shown in figure 4 were obtained to check the convergence of the MSM method to the global minimum. For this reason, some parameters and the tomographic scheme of the model experiment were chosen 'ideal' for tomography problems in order to exclude the influence on the convergence of the MSM method. The number of sources is 24, with six sources located on each side. Wave measurements are taken at each grid point along the perimeter of square $L$ and there are about 100 measurements per wave period in time. The use of broadband sources and receivers with a wide radiation pattern placed at a large number of positions on all sides of the object provide an almost complete set of routes between sources and receivers for sounding the object. This results in high quality image reconstruction.

In addition, all model calculations assumed precisely 'measured' experimental data, i.e. there were no random or systematic errors in the data. The number of iterations of the iterative process was chosen to be greater than 100 at each stage of the MSM method, until the discrepancy stopped decreasing. Of course, such a number of iterations and quality of reconstruction are impossible in processing real experimental data obtained with some error.

Figure 6 shows some results of velocity reconstruction to check the convergence to the global minimum of the MSM method for more realistic values of the experimental parameters. In each experiment in figure 6, a pulse with a central wavelength of 20 mm was used at the first stage of the MSM method, and 5 mm at the second stage. In figure 6(a) the number of sources remains 24, but the step between the measurement points is about 1.5λ ≈ 8 mm and there are about 3–4 measurements per wave period in time. This caused ripples to appear in the image. In figure 6(b), the positions of the measurement points remain as in figure 4, however, the number of sources is reduced to four, with one source located in the center of each side of the square. As can be seen from figure 6(b), this leads to the appearance of artifacts in the image. In figure 6(c), six sources were used, located only on the upper side of the square, but the receivers, as before, were located on all sides of the square. This arrangement of sources results in poor reconstruction at the bottom of the image.

Zoom In Zoom Out Reset image size

In each experiment in figure 6, the quality of the absorption reconstruction is much worse, so it is not shown. As a rule, to obtain the absorption factor with an acceptable quality, it is necessary to increase the number of iterations by several times compared to the number of iterations required for velocity reconstruction. The reconstruction results for other model and real experiments were presented in [10, 39, 44].

The article proves the uniqueness theorem for the solution to the coefficient inverse problem for the wave equation in ${{\boldsymbol{R}}^2}$ with two unknown coefficients: speed of sound and absorption factor.

The original non-linear coefficient inverse problem is reduced to an equivalent system of two uniquely solvable linear integral equations of the first kind with respect to coefficients $c\left( x \right)$ and $a\left( x \right)$. However, the possibility of reducing the original non-linear problem to a linear one depends on the presence of frequencies close to 0 in the spectrum of the sounding signal. Such frequencies are difficult to obtain in physical experiments. In addition, the solution of the equivalent integral equations of the first kind is a highly unstable problem due to the form of the kernel of the integral operator.

Estimations were made to substantiate the MSM for two unknown coefficients. This justification is based on the idea that for small inhomogeneities at sufficiently low frequencies the non-linear inverse problem approaches a linear problem and the residual functional approaches a convex one. The MSM for solving non-linear coefficient inverse problems is not related to the reduction to a linear problem using a limit approach as $w \to 0$, but assumes solving the inverse problem at sufficiently low but non-zero frequencies $w$ instead. For small inhomogeneities which are typical, for example, in medical imaging, conducting experiments at such frequencies does not pose significant problems.

The capabilities of the MSM method are demonstrated on a model inverse problem of determining unknown speed of sound and absorption coefficients in some region. The model problem is focused on ultrasound diagnostics of soft tissues in medicine. The presented model problem demonstrates that in ultrasonic tomography the residual functional can have local minima, and therefore the convergence of gradient iterative minimization methods is not guaranteed. The MSM method effectively solved the problem of non-linearity in two stages.

The paper was published with the financial support of the Ministry of Education and Science of the Russian Federation as part of the program of the Moscow Center for Fundamental and Applied Mathematics under the Agreement № 075-15-2022-284. The research is carried out using the equipment of the shared research facilities of HPC computing resources at Lomonosov Moscow State University.

All data that support the findings of this study are included within the article (and any supplementary files).