Iterative solutions of the inverse problems of frequency sounding and electrical prospecting of layered media

Abstract

The paper presents the iterative solutions of two coefficient inverse problems (CIPs) arising in frequency sounding and electrical prospecting. An iterative algorithm is constructed to obtain such solutions. Exploiting the Beilina–Klibanov approach to CIPs, this algorithm possesses the new iterative and refinement procedures. These features enhance significantly both the spatial and contrast resolutions of reconstructed coefficients. The computational effectiveness of the proposed numerical technique is demonstrated in computational experiments with two applied CIPs: electromagnetic or acoustic frequency sounding and electrical prospecting of layered media. The Slichter–Langer–Tikhonov formulation is exploited as a mathematical model of the latter.

Keywords:

• coefficient inverse problem
• iterative algorithm
• boundary and internal data

AMS Subject Classifications:

• 65J15
• 65L15
• 65N21

1 Introduction

The paper addresses a challenging problem of constructing computationally efficient algorithms for coefficient inverse problems (CIPs), i.e. the inverse problems for the partial differential equations where one needs to determine one or several coefficients from a given solution (or its functionals) on some manifolds (see, e.g. [1–4] and references cited there). In this paper, we present a reconstruction algorithm that produces a sequence of iterates that converges to a good estimate of the sought solution from an initial approximation, which is not necessarily close to this coefficient.

Due to both the non-linearity and ill-posedness of CIPs, providing such a convergence is very important for many, if not for all, reconstruction algorithms. Being applied to a CIP, the traditional non-linear least squares approach normally results in a non-convex minimization problem. As a result, the traditional numerical techniques, such as the gradient or Newton-like methods, may not be applicable or they may fail to converge to the true solution or even to its estimates. On the other hand, the methods of global optimization are extremely time-consuming and heuristic. Therefore, the development of such algorithms is of particular interest to quantitative imaging in many applied areas, such as remote sensing, geophysics, medical diagnostics, non-destructive testing and evaluation, etc.

The goal of this paper is twofold. First, this is to present a modification of the Beilina–Klibanov iterative algorithm [1, 5, 6] for solving some CIPs. Second, this is to demonstrate its computational effectiveness by applying it to frequency sounding and resistivity prospecting of layered media. Introduced by Tikhonov [7] and Cagniard [8] in exploration geophysics, the method of frequency sounding became one of the most popular techniques used in a variety of applications (see, e.g. survey [9]). The term ‘resistivity prospecting’ is associated with geophysical prospecting techniques utilizing measurements of the voltage potential produced by currents injected from the Earth’s surface into the ground (see, e.g. [10]). Unlike the Calderon formulation, which is widely used in electrical impedance tomography with multiple boundary measurements, we utilize the Slichter–Langer–Tikhonov mathematical model [11–13] with a single boundary measurement. Recall that in multiple boundary measurements the number of independent variables in the data is greater than the number of independent variables in the unknown coefficient, whereas in a single boundary measurement they are equal.

To make the underlying idea and technique more descriptive, the one-dimensional (1D) distribution of an unknown coefficient, which is typical for layered media, is considered. However, the algorithm can be generalized to higher dimensions. On the other hand, 1D quantitative imaging itself is widely used in many applications, such as studying Earth interior [14] and seafloor resources,[15] as well as in studying macromolecular biological interactions,[16] nuclear magnetic resonance 1D imaging, mine countermeasures,[6] etc.

One approach to constructing the computational algorithms for CIPs was recently proposed in [2]. It is based on Carleman estimates allowing for constructing some strictly convex least squares functionals at each step of a layer-stripping procedure with respect to the spatial variable. In [1, 6] the layer stripping procedure with respect to the real-valued parameter of the Laplace transform combined with an iterative procedure for treatment of the so-called tails was developed and applied to solving $n$D ($n=1,2,3$) CIPs. Unlike the numerical techniques described in [1, 6], we do not solve an overdetermined boundary value problem when determining the initial approximations. Also, we do not use the marching scheme for performing iterates with respect to the parameter of the Laplace transform. Instead of the WKB (i.e. short wave) approximations, we use the exact solutions for a homogeneous medium in order to initialize the iterative procedure. The main distinctive feature of the proposed algorithm is the refinement procedure, in which we search for a better approximation performing matching the interior field generated by the iterative solution by the solutions of the forward problem in the region of interest. It follows from the results of numerical experiments that these novelties are advantageous. They allow for more accurate estimates of the unknown coefficient.

The paper is formatted as follows. In Section 2, we formulate an inverse problem of frequency sounding. In Section 3, we describe the proposed algorithm. In Section 4, the refinement procedure is outlined and dicussed. In Section 5, we conduct the numerical study of the proposed algorithm by computer simulations of both the frequency sounding and electrical prospecting of layered media. We conclude our investigation in Section 6.

2 An inverse problem of frequency sounding

As in [17], we consider a pulse wave emitted at $z=0$ and propagated through the Epstein layer, which is the half-space $z>0$ filled with the lossless inhomogeneous medium, such that the wave propagation speed is variable for $z\in (0,L)$ and it is constant for $z\ge L$. In this case, the initial boundary value problem for the hyperbolic equation governing the wave propagation is given by1 $$\begin{array}{cc}& c^{-2}(z)(z)U_{tt}-U_{zz}=0z>0,t>0,\hfill \end{array}$$1 2 $$\begin{array}{cc}& U(z,0)=U_t(z,0)=0U(0,t)=\mathit{\delta }(t).\hfill \end{array}$$2 Here, $c(z)$ is the speed of propagation of electromagnetic or acoustic waves. We emphasize that we do not impose the smallness assumption on the coefficient $c(z)$.

Inverse Problem.Let the function$U(z,t)$satisfies the initial boundary value problem (11 $$\begin{array}{cc}& c^{-2}(z)(z)U_{tt}-U_{zz}=0z>0,t>0,\hfill \end{array}$$1 )–(22 $$\begin{array}{cc}& U(z,0)=U_t(z,0)=0U(0,t)=\mathit{\delta }(t).\hfill \end{array}$$2 ). Given the function$U_z(0,t)=f(t),t>0$and parameters$c_0>0,L>0$, such that$c(z)=c_0$for$z\ge L$, find the variable speed$c(z)$on$(0,L)$.

The traditional approach to solving this problem consists of reducing the wave Equation(11 $$\begin{array}{cc}& c^{-2}(z)(z)U_{tt}-U_{zz}=0z>0,t>0,\hfill \end{array}$$1 ) to the Helmholtz equation via the Fourier transform and constructing the globally convergent algorithms. To the author’s knowledge, there exist three such algorithms. The first algorithm (see [18]) is based on the so-called trace (asymptotic) formulae, whereas the second one (see [19]) utilizes the non-linear Riesz transform allowing for reducing the equation to the Volterra integral equation. The third algorithm is based on the convexification method.[2] Unlike these algorithms, we apply the Laplace transform to the problem (11 $$\begin{array}{cc}& c^{-2}(z)(z)U_{tt}-U_{zz}=0z>0,t>0,\hfill \end{array}$$1 ) and(22 $$\begin{array}{cc}& U(z,0)=U_t(z,0)=0U(0,t)=\mathit{\delta }(t).\hfill \end{array}$$2 ). The Laplace transform was also used in [1] for solving the similar inverse problem for the wave equation.

To justify the applicability of the Laplace transform, we assume that $|U(z,t)|\le C·exp({\mathit{\nu }}_{0}t)$ for all $t>0$. Here, $C=const>0$ and ${\mathit{\nu }}_0$ is the given parameter. Then, applying the Laplace transform$$u(x,\mathit{\nu })={\int }_0^{\infty }U(z,t)exp(-\mathit{\nu }t)dt,z>0,\mathit{\nu }>{\mathit{\nu }}_0$$to the problem (11 $$\begin{array}{cc}& c^{-2}(z)(z)U_{tt}-U_{zz}=0z>0,t>0,\hfill \end{array}$$1 ) and (22 $$\begin{array}{cc}& U(z,0)=U_t(z,0)=0U(0,t)=\mathit{\delta }(t).\hfill \end{array}$$2 ) and introducing the dimensionless variables $x=z/L,s=\mathit{\nu }L/c_0,n(x)=c_0/c(Lx)$, one can reduce it to the $s$-parametric family of elliptic boundary value problems3 $$\begin{array}{cc}& u_{xx}(x,s)-s^2{n}^2(x)u(x,s)=0,x\in (0,1),s>s_0,\hfill \end{array}$$3 4 $$\begin{array}{cc}& u(0,s)=1,\hfill \end{array}$$4 5 $$\begin{array}{cc}& u_x(1,s)+su(1,s)=0.\hfill \end{array}$$5 Here, $n(x)=c_0/c(x)$ is the refraction coefficient, where $s_0(s)={\mathit{\nu }}_{0}L/c_0$, and $c(x)$ and $c_0$ are speeds of the wave propagation in the Epstein layer and homogeneous medium. Due to injectivity of the Laplace transform, the existence of a unique solution of (33 $$\begin{array}{cc}& u_{xx}(x,s)-s^2{n}^2(x)u(x,s)=0,x\in (0,1),s>s_0,\hfill \end{array}$$3 )–(55 $$\begin{array}{cc}& u_x(1,s)+su(1,s)=0.\hfill \end{array}$$5 ) follows from classical arguments. Also, for a sufficiently smooth refraction coefficient $u(x,s)>0$ (see, e.g. [2]). In practice, the surface admittance $Y(s)=-u_x(0,s)/u(0,s)$ or impedance is measured in the interval $[\underline{s},\overline{s}],\underline{s}>s_0$. Therefore, introducing by analogy with [2] the new function6 $$\begin{array}{c}\hfill v(x,s)=\frac{1}{s^2}lnu,\end{array}$$6 we obtain the family of problems7 $$\begin{array}{cc}& v_{xx}+s^2{v}_x^2=n^2(x)\text{in}(0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$7 8 $$\begin{array}{cc}& v(0,s)=0,v_x(1,s)=-\frac{1}{s}.\hfill \end{array}$$8 Introducing one more function $q(x,s)={\mathit{\partial }}_{s}v(x,s)$ and differentiating (77 $$\begin{array}{cc}& v_{xx}+s^2{v}_x^2=n^2(x)\text{in}(0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$7 ) and (88 $$\begin{array}{cc}& v(0,s)=0,v_x(1,s)=-\frac{1}{s}.\hfill \end{array}$$8 ) with respect to the $s$-variable, we obtain the mathematical model that does not contain the unknown refraction coefficient $n(x)$9 $$\begin{array}{cc}& q_{xx}+2s^2{q}_x{v}_x=-2s{v}_x^2\text{in}(0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$9 10 $$\begin{array}{cc}& q(0,s)=0,q_x(1,s)=\frac{1}{s^2}.\hfill \end{array}$$10 Note that for every $s\in [\underline{s},\overline{s}]$ the Equation (99 $$\begin{array}{cc}& q_{xx}+2s^2{q}_x{v}_x=-2s{v}_x^2\text{in}(0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$9 ) is integro-differential because of the representation11 $$\begin{array}{c}\hfill v(x,s)=-{\int }_s^{\overline{s}}q(x,\mathit{\mu })d\mathit{\mu }+v(x,\overline{s}).\end{array}$$11 Thus, the inverse problem of frequency sounding is reduced, basically, to estimating the interior field $q(x,s)$ given the data $q_x(0,s)=(2Y(s)/s-Y^{\prime }(s))/s^2\equiv \mathrm{\Psi }(s)$. Indeed, once such an estimate is obtained, the functions $v(x,s)$ and $n(x)$ may also be estimated from (1111 $$\begin{array}{c}\hfill v(x,s)=-{\int }_s^{\overline{s}}q(x,\mathit{\mu })d\mathit{\mu }+v(x,\overline{s}).\end{array}$$11 ) and (77 $$\begin{array}{cc}& v_{xx}+s^2{v}_x^2=n^2(x)\text{in}(0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$7 ).

3 Description of the algorithm

It follows from the previous section that if both functions $q(x,s)$ and $v(x,s)$ are known, then the solution of the inverse problem can be found in the closed form from the Equation (77 $$\begin{array}{cc}& v_{xx}+s^2{v}_x^2=n^2(x)\text{in}(0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$7 ). Another important observation is that if the function $v(x,s)$ is known, then the problem (99 $$\begin{array}{cc}& q_{xx}+2s^2{q}_x{v}_x=-2s{v}_x^2\text{in}(0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$9 ) and (1010 $$\begin{array}{cc}& q(0,s)=0,q_x(1,s)=\frac{1}{s^2}.\hfill \end{array}$$10 ) is linear with respect to $q(x,s)$. We elaborate these observations as follows.

Initialization step   We determine the initial approximations of the functions $v(x,s)$ and $q(x,s)$ as the exact solutions to the boundary value problems (77 $$\begin{array}{cc}& v_{xx}+s^2{v}_x^2=n^2(x)\text{in}(0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$7 )–(1010 $$\begin{array}{cc}& q(0,s)=0,q_x(1,s)=\frac{1}{s^2}.\hfill \end{array}$$10 ) for $n(x)\equiv 1$. That is, we assign$$v^{(0)}(x,s)=-\frac{x}{s},q^{(0)}(x,s)=\frac{x}{s^2}.$$It should be mentioned that the initial approximations determined in [6] from the numerical solution to an auxiliary overdetermined problem are very close to our choice.

Step 1   Beginning with $v^0(x,s)$, we start the iterative process. Assume that the $k$th successive approximation $v^{(k)}(x,s)$ has been determined for all parameters $s\in [\underline{s},\overline{s})$. Then, given $v^{(k)}(x,s)$, we determine $q^k(x,s)$ as follows. For an arbitrary fixed $s\in [\underline{s},\overline{s}]$ denote $L[q](x)=q_{xx}+2s^2{q}_x{v}_x$, where the differential operator from $H^2(0,1)$ to $L_2(0,1)$. Let $F(x)=-2s{v}_x^2$ and consider the Tikhonov functional12 $$\begin{array}{c}\hfill T_{\mathit{\alpha }}[q]={\Vert L[q](x)-F(x)\Vert }_{L_2(0,1)}^2+\mathit{\alpha }{\Vert q\Vert }_{H^2(0,1)}^2,\end{array}$$12 where $\mathit{\alpha }>0$ is the regularization parameter. Since the operator $L[q]$ is linear with respect to $q$, it follows from the theory of ill-posed problems (see, e.g. [20, 21]) that this functional is strictly convex. Then, it follows from the Weierstrass theorem and strict convexity that there exists its unique minimizer $q_{\mathit{\alpha }}(x,s)$ in a set $Q=\{q\in H^2(0,1):q(0,s)=0,q_x(0,s)=\tilde{\mathrm{\Psi }}(s),q_x(1,s)=s^{-2}\}$, and any method of constrained optimization converges to it for every fixed parameter $s\in [\underline{s},\overline{s})$. Then, we assign $q^{k+1}(x,s)=q_{\mathit{\alpha }}(x,s)$. The regularization parameter $\mathit{\alpha }$ is chosen consistently with the level of noise $\mathit{\delta }$ in the ‘measured’ data $\tilde{\mathrm{\Psi }}(s)$, i.e. such that $\Vert \mathrm{\Psi }-\tilde{\mathrm{\Psi }}{\Vert }_{C[\underline{s},\overline{s}]}\le \mathit{\delta }$. In the numerical experiments we used a priori parameter choice $\mathit{\alpha }(\mathit{\delta })={\mathit{\delta }}^2$, where $\mathit{\nu }\in (0,1)$ (see, e.g. [20, 21]).

Step 2   We update the function $v(x,s)$ for every parameter $s\in [\underline{s},\overline{s})$ as follows:13 $$\begin{array}{cc}\hfill v^{(k+1)}(x,s)& =-{\int }_s^{\overline{s}}{q}^{(k+1)}(x,\mathit{\nu })d\mathit{\nu }+v^{(k)}(x,\overline{s}),\hfill \end{array}$$13 14 $$\begin{array}{cc}\hfill v_x^{(k+1)}(x,s)& =-{\int }_s^{\overline{s}}{q}_x^{(k+1)}(x,\mathit{\nu })d\mathit{\nu }+v_x^{(k)}(x,\overline{s}),\hfill \end{array}$$14 15 $$\begin{array}{cc}\hfill v_{xx}^{(k+1)}(x,s)& =-{\int }_s^{\overline{s}}{q}_{xx}^{(k+1)}(x,\mathit{\nu })d\mathit{\nu }+v_{xx}^{(k)}(x,\overline{s}).\hfill \end{array}$$15 Step 3   Let $s^{\ast }\in [\underline{s},\overline{s})$ be a certain number that does not depend on $k$. We update the refraction coefficient as16 $$\begin{array}{c}\hfill n_{k+1}^2(x)=v_{xx}^{(k+1)}\left(x,,,s^{\ast }\right)+s^2{\left(v_x^{(k+1)},\left(x,,,s^{\ast }\right)\right)}^2,\end{array}$$16 where $x\in \left(0,,,1\right)$, and set $n_{k+1}(x)=1$ for $x=0$ and $x=1$.

Step 4   Update the function $v(x,\overline{s})$ by first solving the problem17 $$\begin{array}{cc}& u_{xx}^{(k+1)}(x,\overline{s})-{\overline{s}}^2{n}_{k+1}^2(x)u(x,\overline{s})=0,x\in (0,1),\hfill \end{array}$$17 18 $$\begin{array}{cc}& u^{(k+1)}(0,s)=1,u_x^{(k+1)}(1,s)+s{u}^{(k+1)}(1,s)=0.\hfill \end{array}$$18 and then determining the functions19 $$\begin{array}{cc}& v^{(k+1)}(x,\overline{s})=\frac{ln{u}^{(k+1)}(x,\overline{s})}{{\overline{s}}^2},\hfill \end{array}$$19 20 $$\begin{array}{cc}& v_x^{(k+1)}(x,\overline{s})=\frac{u_x^{(k+1)}(x,\overline{s})}{{\overline{s}}^2{u}^{(k+1)}(x,\overline{s})},\hfill \end{array}$$20 21 $$\begin{array}{cc}& v_{xx}^{(k+1)}(x,\overline{s})=\frac{u_{xx}^{(k+1)}(x,\overline{s})}{{\overline{s}}^2{u}^{(k+1)}(x,\overline{s})}-{\left(\frac{u_x^{(k+1)}(x,\overline{s})}{\overline{s}{u}^{(k+1)}(x,\overline{s})}\right)}^2\hfill \end{array}$$21 for $x\in [0,1]$.

Stopping ruleThe iterative process continues until a stopping rule is fulfilled. If the level of noise $\mathit{\delta }$ in the data $\tilde{\mathrm{\Psi }}(s)$ is given or it can be estimated, then it is natural to stop the process at the iterate22 $$\begin{array}{c}\hfill k_{stop}=min\left\{k:\Vert ,{\mathit{\partial }}_x,u_{n_k},(0,s),-,\tilde{\mathrm{\Psi }},{(s)\Vert }_{L_2[\underline{s},\overline{s}]},\le ,\mathit{\delta }\right\},\end{array}$$22 where the function $u_n(x,s)$ is determined from the solution of the original problem (11 $$\begin{array}{cc}& c^{-2}(z)(z)U_{tt}-U_{zz}=0z>0,t>0,\hfill \end{array}$$1 ) and (22 $$\begin{array}{cc}& U(z,0)=U_t(z,0)=0U(0,t)=\mathit{\delta }(t).\hfill \end{array}$$2 ) with $n(x)=n_k(x)$ for every successive approximation $n_k(x)$. Note that since the norm in (2222 $$\begin{array}{c}\hfill k_{stop}=min\left\{k:\Vert ,{\mathit{\partial }}_x,u_{n_k},(0,s),-,\tilde{\mathrm{\Psi }},{(s)\Vert }_{L_2[\underline{s},\overline{s}]},\le ,\mathit{\delta }\right\},\end{array}$$22 ) represents the residual functional on $n_k(x)$, the number $k_{stop}$ plays the role of the regularization parameter for the entire iterative process. Since in computing we deal with the discrete functions defined on grids, the derivatives are approximated by their finite difference analogous and the integrals in (1313 $$\begin{array}{cc}\hfill v^{(k+1)}(x,s)& =-{\int }_s^{\overline{s}}{q}^{(k+1)}(x,\mathit{\nu })d\mathit{\nu }+v^{(k)}(x,\overline{s}),\hfill \end{array}$$13 )–(1515 $$\begin{array}{cc}\hfill v_{xx}^{(k+1)}(x,s)& =-{\int }_s^{\overline{s}}{q}_{xx}^{(k+1)}(x,\mathit{\nu })d\mathit{\nu }+v_{xx}^{(k)}(x,\overline{s}).\hfill \end{array}$$15 ) – by appropriate quadratures.

Note that if the level $\mathit{\delta }$ is fixed, then any regularizing algorithm, including the proposed one, provides just a certain approximation of the true refraction coefficient $n^{\ast }(x)$. By choosing $\mathit{\alpha }={\mathit{\delta }}^2$ (see, e.g. [20], Chapters 3 and 4), one can easily obtain an error estimate for all $k\in \left[1,,,k_{stop}\right]$23 $$\begin{array}{c}\hfill \Vert n_k(x)-n^{\ast }{(x)\Vert }_{L_2(0,1)}\le {\mathit{\epsilon }}^r,r\in (0,1),\end{array}$$23 where $\mathit{\epsilon }=\mathit{\epsilon }(\mathit{\delta },k_{stop})$ and the number $r$ does not depend on $k,\mathit{\delta },n_k(x)$ and $n^{\ast }(x)$; it does not generally imply the convergence $n_k(x)\to n^{\ast }(x)$ as $k\to \infty$. Therefore, the role of numerical experiments in estimating the approximation error is crucial.

4 Refinement of $n_{k_{stop}}$

It follows from the previous section that the proposed algorithm provides a sufficient proximity of successive approximations to the true coefficient if the level of noise $\mathit{\delta }$ is sufficiently small. In practice, however, it is more typical than exceptional that such a condition is not fulfilled. In this case, the approximation $n_{k_{stop}}(x)$ may not possess high spatial and contrast resolutions, i.e. the quantitative imaging may be problematic. To alleviate this restriction, we propose to take advantage of knowledge of the pair $(\tilde{n}(x),\tilde{u}(x,s))$ everywhere in $(0,1)$, where24 $$\begin{array}{c}\hfill \tilde{u}(x,s)=exp\left[s^2,v^{(k_{stop})},(x,s)\right]\end{array}$$24 and $\tilde{n}(x)=n_{k_{stop}}(x)\in G$, where $G$ is the set of sufficiently smooth functions $n\left(x\right)$ (e.g. $n\in C^2(0,1)$). We observe that although for a certain $s_{\ast }\in [\underline{s},\overline{s})$ the function $\tilde{u}(x,s_{\ast })$ may be close to the Laplace transform of the solution of the original problem (11 $$\begin{array}{cc}& c^{-2}(z)(z)U_{tt}-U_{zz}=0z>0,t>0,\hfill \end{array}$$1 ) and (22 $$\begin{array}{cc}& U(z,0)=U_t(z,0)=0U(0,t)=\mathit{\delta }(t).\hfill \end{array}$$2 ), it does not satisfy it. This observation motivates matching the interior field $\tilde{u}(x,s_{\ast })$ with the Laplace transform $u(x,s_{\ast };n(x))$ of the solution of the original problem (11 $$\begin{array}{cc}& c^{-2}(z)(z)U_{tt}-U_{zz}=0z>0,t>0,\hfill \end{array}$$1 ) and (22 $$\begin{array}{cc}& U(z,0)=U_t(z,0)=0U(0,t)=\mathit{\delta }(t).\hfill \end{array}$$2 ). We are interested in finding a new approximation of $n^{\ast }(x)$ whose error estimate is better than the error estimate for $\tilde{n}(x)$.

To find such an approximation, we consider the injective and continuous map $F:G\to L_2\left(0,,,1\right)$, where $F\left(n\right)=u\left(x,,,s_{\ast },;,n\right),\forall n\left(x\right)\in G$, and minimize the Tikhonov functional25 $$\begin{array}{c}\hfill T_{\mathit{\lambda }}[n]=\Vert F\left(n\right)-\tilde{u}\left(x,,,s_{\ast },;,n\right){\Vert }_{L_2\left(0,,,1\right)}^2+\mathit{\lambda }{\Vert n-\tilde{n}\Vert }_{L_2(0,1)}^2,\mathit{\lambda }>0,\end{array}$$25 where $n_{k_{stop}}\in G$. Because the set $G$ is bounded and compact, the existence of a unique minimizer in a neighbourhood of the function $n^{\ast }(x)$ follows from the Weierstrass theorem and strict convexity of $T_{\mathit{\lambda }}[n]$. Note that the estimate (2323 $$\begin{array}{c}\hfill \Vert n_k(x)-n^{\ast }{(x)\Vert }_{L_2(0,1)}\le {\mathit{\epsilon }}^r,r\in (0,1),\end{array}$$23 ) implies $\Vert u(x,s_{\ast };n^{\ast })-\tilde{u}(x,s_{\ast }){\Vert }_{L_2(0,1)}\le \mathit{\gamma },\mathit{\gamma }>0$.

Let $n_{\mathit{\lambda }}(x)$ be the minimizer of the Tikhonov functional (2525 $$\begin{array}{c}\hfill T_{\mathit{\lambda }}[n]=\Vert F\left(n\right)-\tilde{u}\left(x,,,s_{\ast },;,n\right){\Vert }_{L_2\left(0,,,1\right)}^2+\mathit{\lambda }{\Vert n-\tilde{n}\Vert }_{L_2(0,1)}^2,\mathit{\lambda }>0,\end{array}$$25 ). It follows directly from [22] that if the reqularization parameter is chosen as $\mathit{\lambda }(\mathit{\gamma })={\mathit{\gamma }}^{2\mathit{\nu }},\mathit{\nu }\in (0,1/2)$, then26 $$\begin{array}{c}\hfill ∥,n_{\mathit{\lambda }(\mathit{\gamma })},-,n^{\ast }∥\le \mathit{\xi }\Vert n^{\ast }-n_{k_{stop}}\Vert ,\end{array}$$26 where $\mathit{\xi }\in \left(0,,,1\right)$. This is the motivation of the refinement procedure.

5 Numerical study

In this section, we apply the proposed algorithm to the inverse problems of frequency sounding and electrical prospecting of layered media.

5.1 Computer simulation of frequency sounding

To generate the frequency sounding data $\overline{\mathit{\phi }}(s)$, we solved numerically the problem  (33 $$\begin{array}{cc}& u_{xx}(x,s)-s^2{n}^2(x)u(x,s)=0,x\in (0,1),s>s_0,\hfill \end{array}$$3 )–(55 $$\begin{array}{cc}& u_x(1,s)+su(1,s)=0.\hfill \end{array}$$5 ). To provide a high accuracy of reconstruction, we represented the solution of this problem in the form $u=u_e+\overline{u}$, where $u_e=e^{-sx}$ is its solution for $n(x)\equiv 1$, and considered the problem27 $$\begin{array}{cc}& {\overline{u}}_{xx}(x,s)-s^2{n}^2(x)\overline{u}(x,s)=-s^2{e}^{-sx}(1-n^2(x)),x\in (0,1),\hfill \end{array}$$27 28 $$\begin{array}{cc}& \overline{u}(0,s)=0,\hfill \end{array}$$28 29 $$\begin{array}{cc}& {\overline{u}}_x(1,s)+s\overline{u}(1,s)=0\hfill \end{array}$$29 for every $s\in [\underline{s},\overline{s}]$. The integro-interpolation (or balance) method (see [23], Section 3.2) was used to construct the homogeneous difference scheme for this problem with both the smooth and piecewise-smooth refraction coefficients. Note that if the refraction coefficient is smooth, then we obtain the approximation of order two. The resulted system of linear equations was numerically solved by a stable version of the elimination method (see [23], Section 1.1). Also, since the solution of this difference problem approximates the solution of the problem (2727 $$\begin{array}{cc}& {\overline{u}}_{xx}(x,s)-s^2{n}^2(x)\overline{u}(x,s)=-s^2{e}^{-sx}(1-n^2(x)),x\in (0,1),\hfill \end{array}$$27 )–(2929 $$\begin{array}{cc}& {\overline{u}}_x(1,s)+s\overline{u}(1,s)=0\hfill \end{array}$$29 ), the function $\overline{\mathit{\phi }}(s)$ approximates ${\overline{u}}_x(0,s)$. The level of error of this approximation depends on several factors, such as smoothness of the refraction coefficient, grid, CPU, etc. In the numerical experiments, we utilized the uniform grids containing from 64 to 100 nodes, and the refraction coefficient was assumed to be either smooth or piecewise-constant. Under such conditions, it was estimated from comparison with the model refraction coefficients that the level of noise $\mathit{\delta }$ in the data varied approximately from $5·{10}^{-6}$ (if $n(x)$ is smooth) to $5·{10}^{-4}$ (if $n(x)$ is piecewise-constant). Since $v=\overline{v}-x/s$ and $q=\overline{q}+x/s^2$, where$$\overline{v}=\frac{ln(\overline{u}+e^{-sx})}{s^2}+\frac{x}{s}\text{and}\overline{q}={\mathit{\partial }}_s\overline{v},$$the two-point problem for the function $\overline{q}$ has the form30 $$\begin{array}{cc}& {\overline{q}}_{xx}(x,s)+2s(s{\overline{v}}_x-1){\overline{q}}_x(x,s)=F({\overline{v}}_x;x,s),x\in (0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$30 31 $$\begin{array}{cc}& \overline{q}(0,s)=0,\hfill \end{array}$$31 32 $$\begin{array}{cc}& {\overline{q}}_x(1,s)=0,\hfill \end{array}$$32 where $F({\overline{v}}_x;x,s)=-2{\overline{v}}_x(s{\overline{v}}_x-1)$, and the boundary data is given by33 $$\begin{array}{c}\hfill {\overline{q}}_x(x,0)=\overline{\mathit{\psi }}(s)={\mathit{\partial }}_s\left(\frac{\overline{\mathit{\phi }}(s)}{s^2}\right).\end{array}$$33 In computing, the operator ${\mathit{\partial }}_s$ is understood in the sense of differentiating an interpolating cubic spline of $\overline{\mathit{\phi }}(s)$. The non-linear problem (3030 $$\begin{array}{cc}& {\overline{q}}_{xx}(x,s)+2s(s{\overline{v}}_x-1){\overline{q}}_x(x,s)=F({\overline{v}}_x;x,s),x\in (0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$30 )–(3333 $$\begin{array}{c}\hfill {\overline{q}}_x(x,0)=\overline{\mathit{\psi }}(s)={\mathit{\partial }}_s\left(\frac{\overline{\mathit{\phi }}(s)}{s^2}\right).\end{array}$$33 ) has two unknowns $(\overline{v},\overline{q})$, and it looks similar to the problem (99 $$\begin{array}{cc}& q_{xx}+2s^2{q}_x{v}_x=-2s{v}_x^2\text{in}(0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$9 ) and (1010 $$\begin{array}{cc}& q(0,s)=0,q_x(1,s)=\frac{1}{s^2}.\hfill \end{array}$$10 ). Therefore, we applied the algorithm described in the Section 3. In the iterative process, we updated the refraction coefficient as$$n_{k+1}^2(x)=1+\left(G_k,(x),-,2,s^{\ast },H_k,(x),+,s^{\ast },H_k^2,(x)\right),$$where$$H_k(x)=-{\int }_{s^{\ast }}^{\overline{s}}{\overline{q}}_x^{(k)}(x,\mathit{\nu })d\mathit{\nu }+{\overline{v}}_x^{(k)}(x,\overline{s}),$$$$G_k(x)=-{\int }_{s^{\ast }}^{\overline{s}}{\overline{q}}_{xx}^{(k)}(x,\mathit{\nu })d\mathit{\nu }+{\overline{v}}_{xx}^{(k)}(x,\overline{s}).$$In the numerical experiments, the integrals were approximated by a quadrature formula. Although from a theoretical standpoint, the parameter $s^{\ast }$ can be chosen arbitrarily, in the numerical experiments it was chosen to provide the best accuracy of reconstruction of some benchmark refraction coefficients. The function $\overline{v}(x,\overline{s})$ at $\overline{s}$ was updated as follows. We first solve numerically the two-point problem (2727 $$\begin{array}{cc}& {\overline{u}}_{xx}(x,s)-s^2{n}^2(x)\overline{u}(x,s)=-s^2{e}^{-sx}(1-n^2(x)),x\in (0,1),\hfill \end{array}$$27 )–(2929 $$\begin{array}{cc}& {\overline{u}}_x(1,s)+s\overline{u}(1,s)=0\hfill \end{array}$$29 ) with the $(k+1)$th successive approximation $n_{k+1}(x)$ of the refraction coefficient. Let ${\overline{u}}^{(k+1)}(x,\overline{s})$ be such a solution. Then the updated function $\overline{v}(x,\overline{s})$ and its two derivatives were computed as34 $$\begin{array}{cc}& {\overline{v}}^{(k+1)}(x,\overline{s})=\frac{ln\left({\overline{u}}^{(k+1)}(x,\overline{s})+exp,\left(-,\overline{s},x\right)\right)}{{\overline{s}}^2}+\frac{x}{\overline{s}},\hfill \\ & {\overline{v}}_x^{(k+1)}(x,\overline{s})=\frac{{\overline{u}}_x^{(k+1)}(x,\overline{s})-\overline{s}{e}^{-sx}}{{\overline{s}}^2\left({\overline{u}}^{(k+1)}(x,\overline{s})+exp,\left(-,\overline{s},x\right)\right)}+\frac{1}{\overline{s}},\hfill \\ & {\overline{v}}_{xx}^{(k+1)}(x,\overline{s})=\frac{{\overline{u}}_{xx}^{(k+1)}(x,\overline{s})+{\overline{s}}^2{e}^{-sx}}{{\overline{s}}^2\left({\overline{u}}^{(k+1)}(x,\overline{s})+exp,\left(-,\overline{s},x\right)\right)}-{\left(\frac{{\overline{u}}_x^{(k+1)}(x,\overline{s})-\overline{s}{e}^{-sx}}{\overline{s}\left({\overline{u}}^{(k+1)}(x,\overline{s})+exp,\left(-,\overline{s},x\right)\right)}\right)}^2.\hfill \end{array}$$34 Once the inequality$$\Vert {\mathit{\partial }}_x{\overline{u}}_{n_k}\left(0,,,s\right)-\overline{\mathit{\phi }}{(s)\Vert }_{L_2\left(\underline{s},,,\overline{s}\right)}\le \mathit{\delta }$$was fulfilled, the iterative process was terminated. The approximate solution $n_{k_{stop}}(x)$ was refined as follows. We defined the function35 $$\begin{array}{c}\hfill {\overline{u}}^{(k_{stop})}(x,s_{\ast })=u_0(x,s_{\ast })exp[s_{\ast }^2\overline{v}{}^{(k_{stop})}(x,s_{\ast })]\end{array}$$35 and minimized the functional36 $$\begin{array}{c}\hfill J_{\mathit{\lambda }}(n)=\Vert {\overline{u}}_n\left(x,,,s_{\ast }\right)-{\overline{u}}^{(k_{stop})}(x,s_{\ast }){\Vert }_{L_2(0,1)}^2+\mathit{\lambda }{\Vert n(x)-n_{k_{stop}}(x)\Vert }_{L_2(0,1)}^2.\end{array}$$36 Since the functional $J_{\mathit{\lambda }}(n)$ is strictly convex, the Powell method with the initial approximation ${\overline{n}}_{\mathit{\lambda }}(x)$ was used to find its minimizer.

We first simulated electromagnetic frequency sounding of dry layered soil, which is usually performed by ground penetrating radars for the purpose of land mine detection and classification. We assumed that the relative magnetic permeability and electrical permittivity of the background media are equal to 1, and $n(x)=\sqrt{\mathit{\epsilon }(x)/{\mathit{\epsilon }}_0}$. The realistic values of $n(x)$ (see, e.g. [24]) were used in the numerical experiments. The distribution of the relative electrical permittivity of a subsurface land mine was simulated by a Gaussian curve. Figure 1 shows the typical boundary data $\overline{\mathit{\phi }}(s)$ generated by this mine embedded in the homogeneous background (air) with $\mathit{\epsilon }=1$. In Figure 2, the effect of refinement of ${\mathit{\epsilon }}_{\mathit{\lambda }}(x)$ via fitting the interior data in comparison with refinement via fitting the boundary data is shown. In this experiment, the non-perturbed frequency sounding data were used, and the regularization parameter $\mathit{\lambda }$ providing the best accuracy was chosen $3.6·{10}^{-3}$.

Figure 1. The simulated electromagnetic frequency data $\overline{\mathit{\phi }}(s)$ in the interval $[1,10]$. The noiseless data are shown by bullets and the data corrupted by noise at the level of 5% are shown by asterisks.

Figure 2. Reconstructions of the relative electrical permittivity of a single mid-contrast land mine in the air: (1) without refinement (dashed), (2) with refinement using the boundary data (asterisks) and (3) with refinement using the interior data (bullets). The original distribution is shown by a solid line.

Figure 3. Recovering the relative electrical permittivity of a single mid-contrast land mine from the perturbed data (the level of noise is 5%) without refinement (asterisks) and with refinement (bullets).

Figure 4. Reconstructions of the relative electrical permittivity of a single low-contrast (left), mid-conrast (middle) and high-contrast (right) land mine in the air. The results of reconstruction are shown by asterisks (without refinement) or by bullets (with refinement). The original distribution is shown by a solid line.

Figure 5. Reconstructions of the relative electrical permittivity of a low-contrast land mine embedded in clutter at $x=0.2$ without (asterisks) and with (bullets) refinement. The original distribution is shown by a solid line.

Figure 6. Results of recovering the sound speed profile from the data obtained on the see surface: reconstructions without (asterisks) and with (bullets) refinement. The original sound speed profile is shown by the solid line.

To demonstrate the robustness of the proposed algorithm, we simulated the perturbed data by adding the normally distributed noise to $\overline{\mathit{\phi }}(s)$, i.e. $\tilde{\mathit{\phi }}(s)=\overline{\mathit{\phi }}(s)+\mathit{\gamma }R$, where $R$ is the normally distributed pseudo-random vector with mean zero and standard deviation of 1, and $\mathit{\gamma }$ is the model standard deviation chosen as $\mathit{\gamma }=\mathit{ϵ}\Vert \overline{\mathit{\phi }}\Vert /\Vert R\Vert$, where $\mathit{ϵ}=\Vert \tilde{\mathit{\phi }}-\overline{\mathit{\phi }}\Vert /\Vert \overline{\mathit{\phi }}\Vert$. The perturbed data are shown in Figure 1. Figure 3 shows the mean relative electrical permittivity obtained from a sample containing 25 realization of the pseudo-random vector $R$. Since robustness of the proposed algorithm was also confirmed in all other numerical experiemnts, we limit our next demonstrations by the results of recovering the material properties from the noiseless data.

Figure 7. Results of recovering the sound speed profile from the data obtained on the see floor: reconstructions without (asterisks) and with (bullets) refinement. The original sound speed profile is shown by the solid line.

Since the land mines may be filled with the different explosives and they may be covered by either metallic or plastic materials, their averaged relative electrical permittivity may range from 2 to 25 or more units. Figure 4 shows the results of reconstruction for the different contrasts. In all such cases, the effect of refinement was significant.

In practice, even a near subsurface, where land mines are usually embedded in, is not homogeneous. This may be due to variations of sand and clay porosity, water saturation, presence of abris, etc. These inhomogeneities lead to unwanted variations in GPR signals that are referred to the clutter. We simulated the model clutter containing a continuously inhomogeneous layer with the relative electrical permittivity ranging from 2.4 to 4.75. Such a variation is typical, for instance, for a mixture of sand and clay with the variable porosity and water saturation. We investigated how the presence of clutter affects reconstruction of the relative electrical permittivity of land mines. The reconstructed relative permittivities are shown in Figure 5. Clearly, the presence of clutter affects significantly the quality of quantitative imaging. The higher the contrast the lower the quality, though a near surface land mine may still be detected and classified.

Also, we simulated acoustic frequency sounding of layered marine environments. In these cases, $n(z)=c_0/c(z)$, where $c(z)$ is the sound speed profile, and $c_0$ is the reference sound speed. In the numerical experiments, we used the values of parameters typical for the shallow waters and water-saturated sediments.(see, e.g. [25]) To simulate the sound speed profile in the water column, we used the simplified Wilson’s formula [26]$$c(z)=1449+4.6T(z)-0.055T{(z)}^2+0.0003T{(z)}^3,$$where $T(z)=20-\mathit{\beta }z$ is the temperature (Celcius), $z$ is the depth (m) and $S(z)=34+\mathit{\tau }z$ is the salinity (psu). The sound speed profile in the sea floor sediments was simulated by a piecewise-constant average compressional wave speed in the range from 1685 to 1710 ${\mathrm{m}\mathrm{s}}^{-1}$. Inversion was performed for the dimensionless parameter $s$ ranging from 1 to 10 on a grid containing 64 nodes. Figure 6 shows the results of inversion from the surface of the water column. Clearly, the proposed algorithm performed nicely for the water column, which is of particular interest to underwater acoustics. However, the fine structure of the sediment layer was not resolved. The results of inversion from the sea floor shown in Figure 7 demonstrate the effectiveness of refinement, though the accuracy of quantitative imaging diminished as we advanced into the sediment layer.

5.2 Computer simulation of electrical prospecting

Recovering the subsurface conductivity in a horizontally stratified Earth from the observed surface potential was originated by Slichter [12] and Langer [11]. In [13], Tikhonov established the uniqueness result for their formulation. In [27], the SLT inverse model, as well as the logarithmic derivative of the conductivity function and Hankel transform was used to justify the local strict convexity and smoothness of the residual functional. In this section, we apply the proposed algorithm for the SLT inverse problem.

Assume that the upper-half space is filled with a conductive medium, such that $\mathit{\sigma }=\mathit{\sigma }(z)\ge const.>0$, and there exists a thin layer $(-\mathit{\epsilon },0),\mathit{\epsilon }>0$, such that its conductivity is constant and equals to $\mathit{\sigma }(0)=1$. The lower half-space $z<-\mathit{\epsilon }$ is supposed to be filled with a dielectric. Introduce the cylindrical coordinates $(r,z)$ and consider a point-like current electrode placed at the point $(r_0,z_0):=(0,z_0),z_0\in (-\mathit{\epsilon },0)$. Because of the cylindrical symmetry, the voltage potential $V(r,z)$ generated by this current satisfies the following conditions.37 $$\begin{array}{cc}& r^{-1}{(r{V}_r)}_r+{\mathit{\sigma }}^{-1}(z){(\mathit{\sigma }(z)V_z)}_z=-\mathit{\delta }\left(r\right)\mathit{\delta }(z-z_0),r>0,z>-\mathit{\epsilon },\hfill \end{array}$$37 38 $$\begin{array}{cc}& \mathit{\sigma }(0)V_z(r,-\mathit{\epsilon })=0,\hfill \end{array}$$38 39 $$\begin{array}{cc}& \underset{|(r,z)|\to \infty }{lim}V(r,z)=0.\hfill \end{array}$$39 We formulate the following inverse problem.

Let the function $V(r,z)$ satisfies the boundary value problem (3636 $$\begin{array}{c}\hfill J_{\mathit{\lambda }}(n)=\Vert {\overline{u}}_n\left(x,,,s_{\ast }\right)-{\overline{u}}^{(k_{stop})}(x,s_{\ast }){\Vert }_{L_2(0,1)}^2+\mathit{\lambda }{\Vert n(x)-n_{k_{stop}}(x)\Vert }_{L_2(0,1)}^2.\end{array}$$36 )–(3838 $$\begin{array}{cc}& \mathit{\sigma }(0)V_z(r,-\mathit{\epsilon })=0,\hfill \end{array}$$38 ) and it is associated with the conductivity $\mathit{\sigma }(z)$. However, the trace $V(r,0)$ of $V(r,z)$ on the surface $z=0$ is unknown. Instead, its approximation $\mathrm{\Psi }(r)$, such that $\Vert V(r,0)-\mathrm{\Psi }(r)\Vert \le \mathit{\delta }$, is known. Given $(\mathrm{\Psi }(r),\mathit{\delta },{\mathit{\sigma }}_0,L)$, find an approximation of $\mathit{\sigma }(z)$ in the finite layer $(0,L)$.

Figure 8. The results of inversion of the surface voltage potential with and without refinement are shown by bullets. The original distribution of conductivity for the four-layer medium is shown by the solid line.

Following,[13] we first apply the zero-order Hankel (Fourier-Bessel) transform$$\tilde{V}(z,\mathit{\zeta })={\int }_0^{\infty }V(r,z)J_0(\mathit{\zeta }r)rdr,$$to the problem (3636 $$\begin{array}{c}\hfill J_{\mathit{\lambda }}(n)=\Vert {\overline{u}}_n\left(x,,,s_{\ast }\right)-{\overline{u}}^{(k_{stop})}(x,s_{\ast }){\Vert }_{L_2(0,1)}^2+\mathit{\lambda }{\Vert n(x)-n_{k_{stop}}(x)\Vert }_{L_2(0,1)}^2.\end{array}$$36 )–(3838 $$\begin{array}{cc}& \mathit{\sigma }(0)V_z(r,-\mathit{\epsilon })=0,\hfill \end{array}$$38 ) and introduce the new variables40 $$\begin{array}{c}\hfill t(z)={\int }_{z_0}^z{\mathit{\sigma }}^{-1}(\mathit{\xi })d\mathit{\xi },\end{array}$$40 $x=t/T,s=\mathit{\sigma }(L)\mathit{\zeta }T,n(x)=\mathit{\sigma }(z(xT))/\mathit{\sigma }(L)$, where $a(t)=\mathit{\sigma }(z(t))$ and$$T={\int }_{z_0}^L{\mathit{\sigma }}^{-1}(\mathit{\xi })d\mathit{\xi }.$$Then, the problem (3636 $$\begin{array}{c}\hfill J_{\mathit{\lambda }}(n)=\Vert {\overline{u}}_n\left(x,,,s_{\ast }\right)-{\overline{u}}^{(k_{stop})}(x,s_{\ast }){\Vert }_{L_2(0,1)}^2+\mathit{\lambda }{\Vert n(x)-n_{k_{stop}}(x)\Vert }_{L_2(0,1)}^2.\end{array}$$36 )–(3838 $$\begin{array}{cc}& \mathit{\sigma }(0)V_z(r,-\mathit{\epsilon })=0,\hfill \end{array}$$38 ) is reduced to the two-point problem41 $$\begin{array}{cc}& u_{xx}(x,s)-s^2{n}^2(x)u(x,s)=0,x\in (0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$41 42 $$\begin{array}{cc}& u_x(0,s)=-\mathit{\delta }(x-x_0),\hfill \end{array}$$42 43 $$\begin{array}{cc}& u_x(1,s)+su(1,s)=0.\hfill \end{array}$$43 Here, $u(x,s)=\tilde{V}(z(xT),s/\mathit{\sigma }(L)T)$. We observe the similarity of this problem with the problem (33 $$\begin{array}{cc}& u_{xx}(x,s)-s^2{n}^2(x)u(x,s)=0,x\in (0,1),s>s_0,\hfill \end{array}$$3 )–(55 $$\begin{array}{cc}& u_x(1,s)+su(1,s)=0.\hfill \end{array}$$5 ) modelling frequency sounding of layered media. The only difference is the surface boundary conditions (4141 $$\begin{array}{cc}& u_{xx}(x,s)-s^2{n}^2(x)u(x,s)=0,x\in (0,1),s\in [\underline{s},\overline{s}],\hfill \end{array}$$41 ), which is due to the presence of the current electrode. This observation motivates the following reformulation of the inverse problem.

Let the function$u(x,s)$be a solution of the problem(4040 $$\begin{array}{c}\hfill t(z)={\int }_{z_0}^z{\mathit{\sigma }}^{-1}(\mathit{\xi })d\mathit{\xi },\end{array}$$40 )–(4242 $$\begin{array}{cc}& u_x(0,s)=-\mathit{\delta }(x-x_0),\hfill \end{array}$$42 ) that corresponds to$n(x)$, but$u(0,s)$is unknown. Instead, an approximation$\mathit{\phi }(s)$, such that$\Vert u(0,s)-\mathit{\phi }(s)\Vert \le {\mathit{\delta }}_{\mathit{\phi }}$, where${\mathit{\delta }}_{\mathit{\phi }}>0$and$\mathit{\phi }(s)=\mathit{\psi }(s/\mathit{\sigma }(L)T)$$$\mathit{\psi }(\mathit{\tau })={\int }_0^{\infty }\mathrm{\Psi }(r)J_0(\mathit{\tau }r)rdr$$is given, find an approximation of$n(x)$in$(0,1)$.

Since the SLT model (4040 $$\begin{array}{c}\hfill t(z)={\int }_{z_0}^z{\mathit{\sigma }}^{-1}(\mathit{\xi })d\mathit{\xi },\end{array}$$40 )–(4242 $$\begin{array}{cc}& u_x(0,s)=-\mathit{\delta }(x-x_0),\hfill \end{array}$$42 ) is similar to the frequency sounding model (33 $$\begin{array}{cc}& u_{xx}(x,s)-s^2{n}^2(x)u(x,s)=0,x\in (0,1),s>s_0,\hfill \end{array}$$3 )–(55 $$\begin{array}{cc}& u_x(1,s)+su(1,s)=0.\hfill \end{array}$$5 ), the specific procedure for reconstructing the conductivity is the same as described in the Section 5.1. However, since the travel time transformation (3939 $$\begin{array}{cc}& \underset{|(r,z)|\to \infty }{lim}V(r,z)=0.\hfill \end{array}$$39 ) maps a uniform grid ${\{z_i\}}_{i=0}^n$ onto a non-uniform grid ${\{t_i\}}_{i=0}^n$, such that$$t_i={\mathit{\sigma }}^{-1}(z_i)h+t_{i-1},t_0=0,t_n=T,$$one needs to transform the reconstructed refraction coefficient back to the conductivity distribution. Let ${\tilde{n}}_i$ be the values of the reconstructed refraction coefficient on an arbitrary uniform grid ${\{{\widehat{x}}_i\}}_{i=0}^n$. Then, to obtain the values ${\tilde{\mathit{\sigma }}}_i$ of the reconstructed conductivity on ${\{z_i\}}_{i=0}^n$, we interpolated the values ${\tilde{n}}_i\mathit{\sigma }(L)$ from the uniform grid ${\{{\widehat{x}}_i\}}_{i=0}^n$ onto the non-uniform grid ${\{t_i/T\}}_{i=0}^n$ corresponding to the uniform grid ${\{z_i\}}_{i=0}^n$. In this case, the ‘optimal’ dimensionless range of the parameter $s$, i.e. $[\underline{s},\overline{s}]$, was chosen as $[0.01,0.1]$. It is demonstrated in Figure 8 that the refinement procedure does not provide the improvement of ${\mathit{\sigma }}_{k_{stop}}(z)$. This is due to the diffusive character of the stationary electric field described by the SLT model (4040 $$\begin{array}{c}\hfill t(z)={\int }_{z_0}^z{\mathit{\sigma }}^{-1}(\mathit{\xi })d\mathit{\xi },\end{array}$$40 )–(4242 $$\begin{array}{cc}& u_x(0,s)=-\mathit{\delta }(x-x_0),\hfill \end{array}$$42 ).

6 Conclusions

We modified the Beilina–Klibanov iterative algorithm for quantitative imaging of layered media. Compared to the original algorithm, it possesses the new iterative and refinement procedures. The modified algorithm was tested in several numerical experiments with the data simulated for electromagnetic and acoustic frequency sounding, as well as for electrical prospecting of layered media. The results of testing demonstrated the computational effectiveness of the modified algorithm when recovering the permittivity of the near surface targets from the electromagnetic frequency sounding data. However, in cases of acoustic sounding and electrical prospecting of geophysical layered configurations, the modified algorithm exhibited the so-called smoothing property (see Figures 6–8). Apparently, this is due to the use of the $L_2$-norm in the Tikhonov functional (2525 $$\begin{array}{c}\hfill T_{\mathit{\lambda }}[n]=\Vert F\left(n\right)-\tilde{u}\left(x,,,s_{\ast },;,n\right){\Vert }_{L_2\left(0,,,1\right)}^2+\mathit{\lambda }{\Vert n-\tilde{n}\Vert }_{L_2(0,1)}^2,\mathit{\lambda }>0,\end{array}$$25 ). Therefore, applying the total variation regularization to these problems will be the subject in our further study.