Analysis of new direct sampling indicators for far-field measurements

TL;DR

This paper analyzes three direct sampling indicators for reconstructing scatterers from far-field acoustic data, demonstrating their equivalence and introducing a new indicator based on Tikhonov regularization, with numerical comparisons.

Contribution

It introduces a novel indicator based on Tikhonov regularization and establishes the equivalence of existing indicators for far-field scattering reconstruction.

Findings

The new Tikhonov-based indicator performs comparably to existing methods.

The indicators studied are mathematically equivalent under certain conditions.

Numerical examples illustrate the effectiveness of the proposed approach.

Abstract

This article focuses on the analysis of three direct sampling indicators which can be used for recovering scatterers from the far-field pattern of time-harmonic acoustic measurements. These methods fall under the category of sampling methods where an indicator function is constructed using the far-field operator. Motivated by some recent work, we study the standard indicator using the far-field operator and two indicators derived from the factorization method. We show the equivalence of two indicators previously studied as well as propose a new indicator based on the Tikhonov regularization applied to the far-field equation for the factorization method. Finally, we give some numerical examples to show how the reconstructions compare to other direct sampling methods.

Figures

Equations

$$z⟼\left|(F{\varphi }_z,{\varphi }_z{)}_{L^2(\mathbb{S})}\right|$$

$$\Delta u^s+k^2{u}^s=0\text{ in }{\mathbb{R}}^d\setminus \overline{D}$$

$$u=u^s+u^i$$

$$\underset{r\to \infty }{\lim }{r}^{(d-1)/2}\left({\partial }_r{u}^s-\mathrm{i}k{u}^s\right)=0.$$

$$u^s(x,\hat{y})=\gamma \frac{{\mathrm{e}}^{\mathrm{i}k|x|}}{|x{|}^{(d-1)/2}}\left\{u^{\infty }(\hat{x},\hat{y})+\mathcal{O}\left(\frac{1}{|x|}\right)\right\}\text{ as }|x|\to \infty$$

$$\gamma =\frac{{\mathrm{e}}^{\mathrm{i}\pi /4}}{\sqrt{8\pi k}}\text{in}{\mathbb{R}}^2\text{and}\gamma =\frac{1}{4\pi }\text{in}{\mathbb{R}}^3$$

$$(Fg)(\hat{x})={\int }_{\mathbb{S}}{u}^{\infty }(\hat{x},\hat{y})g(\hat{y})\mathrm{d}s(\hat{y}),\text{ where }\mathbb{S}=\text{unit sphere/circle.}$$

$$\left(|F|^{1/2}\phi_{z},\phi_{z}\right)_{L^{2}(\mathbb{S})}\quad\text{and}\quad\big{\|}{\color[rgb]{0,0,0}P_{\alpha,\varepsilon}}(|F|)\phi_{z}\big{\|}^{2}_{L^{2}(\mathbb{S})}\,,\quad\text{where}\quad\phi_{z}=\mathrm{e}^{-\mathrm{i}k\hat{x}\cdot z}\,.$$

$$|F{|}^{1/2}{g}_z={\varphi }_z\text{ for }z\in {\mathbb{R}}^d,$$

$$H:L^{2}(\mathbb{S})\to\big{[}L^{2}(D)\big{]}^{d+1}\quad\text{ is given by }\quad Hg=\left(\nabla{v}_{g},{v}_{g}\right)^{\top}$$

$$v_g(x)={\int }_{\mathbb{S}}{\mathrm{e}}^{\mathrm{i}kx\cdot \hat{y}}g(\hat{y})\mathrm{d}s(\hat{y}).$$

$$(v_{\phi_{z}})(x)=\int_{\mathbb{S}}\mathrm{e}^{-\mathrm{i}k{(z-x)}\cdot\hat{y}}\,\mathrm{d}s(\hat{y})=\left\{\begin{array}[]{lr}2\pi J_{0}(k|x-z|)&\,\text{in}\,\,\mathbb{R}^{2}\,,\\ &\\ 4\pi j_{0}(k|x-z|)&\,\text{in}\,\,\mathbb{R}^{3}\,,\end{array}\right.$$

$$\left|(F{\varphi }_z,{\varphi }_z{)}_{L^2(\mathbb{S})}\right|$$

$$\left|(F\phi_{z},\phi_{z})_{L^{2}(\mathbb{S})}\right|=\mathcal{O}\left(\mathrm{dist}(z,D)^{1-d}\right){\color[rgb]{0,0,0}\quad\text{ as }\quad\mathrm{dist}(z,D)\to\infty.}$$

$$W_{\text{FDSM}}(z)={\left(|F{|}^{1/2}{\varphi }_z,{\varphi }_z\right)}_{L^2(\mathbb{S})}\text{ and }{W}_{\text{DSM}}(z)=\left|(F{\varphi }_z,{\varphi }_z{)}_{L^2(\mathbb{S})}\right|$$

$$S=I+2\mathrm{i}k|\gamma {|}^{2}F.$$

$$Fg=\sum _{j=1}^{\infty }{\lambda }_j(g,{\psi }_j{)}_{{}_{L^2(\mathbb{S})}}{\psi }_j\text{for all }g\in L^2(\mathbb{S}).$$

$$F=|F{|}^{1/2}Q|F{|}^{1/2}$$

$$|F{|}^{1/2}g=\sum _{j=1}^{\infty }\sqrt{|{\lambda }_j|}(g,{\psi }_j{)}_{{}_{L^2(\mathbb{S})}}{\psi }_j\text{ and }Qg=\sum _{j=1}^{\infty }\frac{{\lambda }_j}{|{\lambda }_j|}(g,{\psi }_j{)}_{{}_{L^2(\mathbb{S})}}{\psi }_j.$$

$$\displaystyle\big{|}(Qg,g)_{{L^{2}(\mathbb{S})}}\big{|}\geq\mu\|g\|^{2}_{{L^{2}(\mathbb{S})}}\quad\text{ and }\quad\|Qg\|_{{L^{2}(\mathbb{S})}}\leq\|g\|_{{L^{2}(\mathbb{S})}}$$

$$\left|(F{\varphi }_z,{\varphi }_z{)}_{L^2(\mathbb{S})}\right|$$

$$\Vert {\varphi }_z{\Vert }_{L^2(\mathbb{S})}^2=2^{d-1}\pi$$

$$\Vert Tv{\Vert }_{\mathcal{V}}^2\le \Vert T\Vert (Tv,v{)}_{\mathcal{V}}\text{ for all }v\in \mathcal{V}.$$

$$\left|(F{\varphi }_z,{\varphi }_z{)}_{L^2(\mathbb{S})}\right|$$

$$\frac{\mu }{2^{d-1}\pi }{W}_{\text{FDSM}}^2(z)\le W_{\text{DSM}}(z)\le \Vert |F{|}^{1/2}\Vert W_{\text{FDSM}}(z).$$

$$c_1{\left|{\left(|F{|}^{1/2}{\varphi }_z,{\varphi }_z\right)}_{L^2(\mathbb{S})}\right|}^2\le \left|(F{\varphi }_z,{\varphi }_z{)}_{L^2(\mathbb{S})}\right|\le c_2{\left(|F{|}^{1/2}{\varphi }_z,{\varphi }_z\right)}_{L^2(\mathbb{S})}.$$

$$\left|\left(|F|^{1/2}\phi_{z},\phi_{z}\right)_{L^{2}(\mathbb{S})}\right|^{2}=\mathcal{O}\left(\mathrm{dist}(z,D)^{1-d}\right){\color[rgb]{0,0,0}\quad\text{ as }\quad\mathrm{dist}(z,D)\to\infty.}$$

$$W_{\text{FDSM}}(z)=\sum _{j=1}^{\infty }\sqrt{|{\lambda }_j|}{\left|({\varphi }_z,{\psi }_j{)}_{{}_{L^2(\mathbb{S})}}\right|}^2$$

$$z⟼\sum _{j=1}^{\infty }\frac{1}{|{\lambda }_j|}{\left|({\varphi }_z,{\psi }_j{)}_{{}_{L^2(\mathbb{S})}}\right|}^2.$$

$$\Delta u^s+k^2{u}^s=0\text{ in }{\mathbb{R}}^d\setminus \overline{D}\text{ and }\mathcal{B}(u^s)=-\mathcal{B}(u^i)\text{ on }\partial D$$

$$\mathcal{B}(\phi )=\phi \text{ or }\mathcal{B}(\phi )={\partial }_{\nu }\phi +\gamma (x)\phi$$

Full text

Keywords: Inverse scattering, sampling methods, far-field pattern, factorization method, Tikhonov regularization

AMS subject classifications: 35J05, 78A46

1 Introduction

In this paper, the inverse problem of recovering the shape of an unknown obstacle from far-field measurements is considered. Such a problem arises in many physical applications where one needs to detect structures in a given medium. This has many applications such as medical imaging or detecting defects in complex material structures. An important question is if one can construct fast stable algorithms to reconstruct the obstacle with little to no a priori information. One way to achieve this is to employ a qualitative method (otherwise known as non-iterative or direct methods) such as the linear sampling method (LSM) which was first proposed in [10] or the factorization method (FM) which has been applied to acoustic, electromagnetic, and electrostatic problems in [18] as well as the generalized linear sampling method (GLSM) introduced in [2] which combines theoretical elements of both the LSM and FM. These methods have been used to solve multiple inverse shape problems for elliptic [21], hyperbolic [13], and parabolic [14] systems. The LSM studies the far-field equation and uses the fact that the (regularized) solution to the far-field equation should become unbounded in the exterior of the obstacle. The FM gives that the far-field pattern of the fundamental solution to the Helmholtz equation is in the range of a self-adjoint compact operator defined by the far-field operator if and only if the sampling point is in the region of interest. Therefore, appealing to Picard’s Criterion gives an indicator that becomes unbounded in the exterior of the obstacle by using the eigensystem of a known operator. There has also been some work in using numerical regularization techniques when applying the FM in [3]. Recently in [7] and [12] the FM has been extended to reconstructing defects in a known inhomogeneous media. The GLSM studies the minimizer of a functional to solve the far-field equation used in the LSM involving the far-field operator and it is shown that the penalty term in the functional becomes unbounded in the exterior of the obstacle. In the GLSM the penalty term can be defined by using the operator from the FM. All of these methods allow one to construct an indicator function $W(z)$ to recover the unknown obstacle where $W(z)$ is positive in the interior of the obstacle and is (approximately) zero on the exterior of the obstacle. Therefore, to reconstruct the obstacle one can plot the computable function $W(z)$ in a region where the unknown scatterer is assumed to be located. Lately, there has been some interest in analyzing a so-called direct sampling method (DSM) using the theoretical framework used in the FM for far-field data (see for e.g. [22] and [15]). These methods have also been studied for there applicability to reconstruct scatterers from near-field data [9]. The DSM recovers objects by constructing an indicator function by evaluating an inner-product with the measured far-field operator and the known far-field pattern of the fundamental solution to Helmholtz equation. This allows one to prove that these indicators are stable and computationally cheap to implement. Similar to the other sampling methods one wishes to show that the DSM indicator has a specific behavior outside the region of interest, where instead of blowing up outside the obstacle it takes smaller values in the exterior.

Recently, in [22] the DSM using the symmetric factorization of the far-field operator for time-harmonic acoustic data has been studied and the indicator function

[TABLE]

is proposed, where $F$ is the far-field operator and $\phi_{z}$ is the far-field pattern of the fundamental solution located at $z$ to Helmholtz equation. The analysis in this paper works for sound soft/hard, impedance obstacles as well as penetrable isotropic scatterers. In [22] it is shown that the indicator is strictly positive and decays as $|z|\to\infty$. The decay is given by bounding the indicator by the Bessel functions evaluated at $|z-x|$ where $x$ is in the obstacle or on the boundary of the obstacle. It is also shown that the indicator is equivalent to reverse time migration [9] and the orthogonal sampling method [23]. Similar analysis was employed in [16] to validate the DSM for elastic waves as well as the study of its applicability for limited-aperture data. Unlike the other sampling methods the DSM does not require solving an ill-posed problem or minimizing a functional at each sampling point which makes it computationally cheaper than its counterparts. Stability results can be given by simple calculations using basic analytical tools. This makes these methods advantageous to use in order to recover unknown obstacles.

Here we propose the use of two indicator functions based on the FM. This research is motivated by the recent paper [15] where similar indicators are introduced and analyzed for the multi-static response matrix which is the discrete version of the far-field operator. The indicators in [15] are based on the theoretical setting of the FM and one of the proposed indicators is based on a Neumann series approximation to the FM equation. In this paper, we will analyze the indicator function which uses the FM operator and prove that the indicator is equivalent to DSM studied in [22] and derive a new method based on the Tikhonov regularization. The main idea for deriving the new DSM is based on the factorization of the far-field operator given by $F$ which uses $|F|^{1/2}=(F^{*}F)^{1/4}$. If we let the unknown obstacle be denoted by $D$, then the FM gives that $\phi_{z}\in\mathrm{Range}(|F|^{1/2})\iff z\in D$ and our method looks at an approximation of the Tikhonov regularized solution operator of the equation $|F|^{1/2}g_{z}=\phi_{z}\,.$

The rest of the paper is structured as follows. First, we provide a precise problem statement in Section 2 and give some preliminary result. In Section 3 we analyze the first indicator which is defined by an inner-product involving $|F|^{1/2}$ and $\phi_{z}$ which was originally considered in [15]. The authors of [15] did not succeed in showing that the indicator is equivalent to the one proposed in [22] and here we will show that the indicators are equivalent. In Section 4 we expand on the idea in [15] to use an approximate solution operator to define a new indicator function. In order to define our approximate solution operator we consider using Tikhonov regularization. Section 5 is devoted to giving some numerical examples of recovering scatterers in two and three dimensions as well as comparing the new indicators with the DSM indicators to show that our reconstructions can compete and in some cases outperform the DSM indicator considered in [22]. In our experiments we see that the reconstructions take seconds to compute making these methods computationally cheap to implement and analytically rigorous. Lastly, a short summary is given in Section 6.

2 Problem statement

In our analysis, we will consider the time-harmonic acoustic scattering problem for anisotropic media. The LSM and FM have been applied to the inverse scattering problem of recovering the scatterer from the measured far-field pattern in [6] and [19], respectively. We now derive two DSMs for this inverse shape problem. To this end, we formulate the direct scattering problems under consideration in $\mathbb{R}^{d}$ for $d=2$ or $d=3$. The scatterer $D\subset\mathbb{R}^{d}$ may be made up of multiple simply connected components with $C^{2}$ boundary $\partial D$ having unit outward normal $\nu$. We consider the scattering by a plane incident wave $u^{i}(x,\hat{y})=\mathrm{e}^{\mathrm{i}k{x}\cdot\hat{y}}$ for a given incident direction $\hat{y}\in\mathbb{R}^{d}$ such that $|\hat{y}|=1$ and wave number $k>0$. This gives that the radiating scattered field $u^{s}(x,\hat{y})$ and the total field $u=u^{s}+u^{i}\in H^{1}_{loc}(\mathbb{R}^{d})$ satisfy the boundary value problem

[TABLE]

along with the Sommerfeld radiation condition

[TABLE]

Here the normal and conormal derivative on the boundary is given by ${\partial_{\nu}\varphi}=\nu\cdot\nabla\varphi$ and ${\partial_{\nu_{A}}\varphi}=\nu\cdot A\nabla\varphi$, respectively, where $\varphi$ is a sufficiently regular function defined on $\partial D$. The Sommerfeld radiation condition is assumed to be satisfied uniformly with respect to $\hat{x}=x/r$, $r=|x|\,.$

The real-valued coefficient matrix $A(x)\in C^{1}(D,\mathbb{R}^{d\times d})$ is symmetric uniformly positive definite in $D$ and the scalar function $n(x)\in C(D)$ is real-valued, denotes the material parameters of the obstacle $D$ where we assume that $I-A$ and $1-n$ are supported in $D$. Under these assumptions we have that the scattering problem (1)–(2) with the radiation condition is well-posed for all $\hat{y}$. It can be shown that the radiating scattered field $u^{s}$ has the expansion

[TABLE]

where the constant

[TABLE]

with $u^{\infty}(\hat{x},\hat{y})$ being the far-field pattern depending on the incident direction $\hat{y}$ and the observation direction $\hat{x}$. We now define the far-field operator $F:L^{2}(\mathbb{S})\longmapsto L^{2}(\mathbb{S})$

[TABLE]

Since the far-field pattern is analytic it is clear that $F$ is a compact operator. The two indicators we study are given by

[TABLE]

Since $F^{*}F$ is a positive self-adjoint compact operator we can define $|F|^{1/2}=(F^{*}F)^{1/4}$ via the spectral decomposition. Here, $P_{\alpha,\varepsilon}(t)$ is a polynomial defined on the interval $[0,\|F\|]$ that approximates the Tikhonov regularized filter function with regularization parameter ${\alpha}$ and accuracy $\varepsilon$ of the equation

[TABLE]

see Section 4 for details. The merits of developing these new indicators is to theoretically justify the numerical investigation in [15] for the continuous (and discretized) setting. We also wish to analyze the new indicators to connect the recent developments for DSMs to earlier qualitative reconstructive algorithms. In our experiments computing the singular values/vector as well as constructing the polynomial $P_{\alpha,\varepsilon}$ adds a negligible increase in computing time for a modest sized far-field matrix.

Before we begin, we will show that the indicator $\left|(F\phi_{z},\phi_{z})_{L^{2}(\mathbb{S})}\right|$ decays as dist$(z,D)$ increases for an anisotropic scatterer (see [22] for other scatterers). It has been shown in [5] that $F$ has the factorization $F=H^{*}TH$ such that

[TABLE]

where $v_{g}$ is the Herglotz wave function defined as

[TABLE]

Here $H^{*}$ is the adjoint operator to $H$ and $T$ is a bounded linear operator from \big{[}L^{2}(D)\big{]}^{d+1} to itself. We have the identity

[TABLE]

where $J_{0}$ is the zeroth order Bessel function of the first kind and $j_{0}$ the zeroth order spherical Bessel function of the first kind (see for e.g. [22]). The factorization of the far-field operator implies that

[TABLE]

Recall that $J_{0}(t)$ and its derivatives decay at a rate of $t^{-1/2}$ as $t\to\infty$. In $\mathbb{R}^{3}$ we have that $j_{0}(t)$ and its derivatives decay rate of $t^{-1}$ as $t\to\infty$. Therefore, the above inequality and the decay of the Bessel functions gives that

Theorem 2.1**.**

For all $z\in\mathbb{R}^{d}\setminus\overline{D}$

[TABLE]

We then see that the indicator will decay as $z$ moves away from the scatterer $D$. Note that Theorem 2.1 is valid for the case when $A$ and $n$ are complex-valued functions.

3 A factorization based direct sampling method

In this section, we will study the indicator using the operator $|F|^{1/2}=(F^{*}F)^{1/4}$. It is well known that one can uniquely recover the scatterer $D$ using the far-field pattern (see Chapter 6 of [4]). The purpose of analyzing the DSM is to derive stable and computationally simple reconstruction algorithms. The stability of the DSMs proposed in this section have been studied in [22] and [15] (for the discretized case). We will show that the corresponding indicator functions

[TABLE]

with $\phi_{z}=\mathrm{e}^{-\mathrm{i}k\hat{x}\cdot z}$ are equivalent. This crucial theoretical result is needed to prove the validity of the new indicator $W_{\text{{FDSM}}}(z)$ which is not established in [15]. Here $W_{\text{{FDSM}}}(z)$ is the DSM based on the FM and $W_{\text{{DSM}}}(z)$ is the standard DSM studied in [22]. To this end, we will bound $W_{\text{DSM}}(z)$ from above and below by our new indicator $W_{\text{{FDSM}}}(z)$. Since $W_{\text{{DSM}}}(z)$ decays as the dist$(z,D)$ increases this suggests that we can plot either function to recover the scatterer.

To begin, we need a few results for the far-field operator defined in (3) that is associated with (1)–(2) along with the radiation condition. We introduce the scattering operator associated with this direct scattering problem (1)–(2). The scattering operator ${S}:L^{2}(\mathbb{S})\to L^{2}(\mathbb{S})$ is defined by

[TABLE]

Since $A(x)$ and $n(x)$ are assumed to be real-valued, we have that the scattering operator is unitary, ${S}{S}^{*}={S}^{*}{S}=I$ (see for e.g. [4]). Here $I$ denotes the identity operator on $L^{2}(\mathbb{S})$. This implies that the corresponding far-field operator $F:L^{2}(\mathbb{S})\to L^{2}(\mathbb{S})$ is normal and compact. We can now conclude that $F$ has an orthonormal eigenvalue decomposition $(\lambda_{j},\psi_{j})\in\mathbb{C}\times L^{2}(\mathbb{S})$ such that

[TABLE]

Since $F$ is a compact operator we have that $|\lambda_{j}|\to 0$ as $j\to\infty$. Provided that the wave number $k$ is not an interior transmission eigenvalue we have that $F$ is injective with a dense range (see [4, 5]). This implies that $|\lambda_{j}|\neq 0$ for all $j$ and that the set $\{\psi_{j}\}$ is a complete orthonormal set in $L^{2}(\mathbb{S})$. Following the results of [18, Chapter 1] it can be shown that the far-field operator has the factorization

[TABLE]

where we use the eigensystem to define the operators

[TABLE]

It is clear that $|F|^{1/2}$ and $Q$ are bounded linear operators that map $L^{2}(\mathbb{S})$ into itself. Since $F$ is injective we have that $|F|^{1/2}$ is a positive self-adjoint compact operator. Since the scattering operator is unitary we can conclude that for all $g\in L^{2}(\mathbb{S})$

[TABLE]

for some positive constant $\mu$ (see for e.g. Theorem 7.29 in [4] for details).

In [15] the authors did not succeed in bounding $W_{\text{{FDSM}}}(z)$ above by $W_{\text{{DSM}}}(z)$ for the discretized far-field operator which is what is needed for showing the equivalence of these two indicators. From the factorization of the far-field operator (4) and the above inequality (5) we have the following estimates

[TABLE]

since

[TABLE]

holds. The above estimate gives that $W^{2}_{\text{{FDSM}}}(z)$ will have the same decay as dist$(z,D)$ increases which is the critical piece missing in the manuscript [15]. We have just about all we need to prove that the two indicators are equivalent. In order to prove the upper bound we need one last theoretical lemma for positive operators on a Hilbert space to prove the equivalence.

Lemma 3.1**.**

Let $T:\mathcal{V}\to\mathcal{V}$ be a bounded positive operator on a Hilbert space $\mathcal{V}$ then we have that

[TABLE]

Proof.

For the proof we refer the reader to [24, Lemma 2.1]. ∎

Notice that since $F$ is injective we have that $|F|^{1/2}$ is a positive operator acting on the Hilbert space $L^{2}(\mathbb{S})$. Using the above lemma we have that

[TABLE]

Thus, we have proven that $W_{\text{{DSM}}}(z)$ is bounded above by $W_{\text{{FDSM}}}(z)$ which implies that they are equivalent in the sense that

[TABLE]

Therefore, we have that the following result.

Theorem 3.1**.**

There exist two positive constants $c_{1}$ and $c_{2}$ such that

[TABLE]

Since we have shown in Theorem 2.1 that $W_{\text{DSM}}(z)$ is bounded by the squared $H^{1}(D)$ norm of the zeroth order Bessel function in $\mathbb{R}^{d}$ for $d=2$ or $3$ we have that $W_{\text{{DSM}}}(z)=\mathcal{O}\left(\mathrm{dist}(z,D)^{1-d}\right)\,$ as $\,\mathrm{dist}(z,D)\to\infty$ when the sampling point $z\in\mathbb{R}^{d}\setminus\overline{D}$ which gives the following result.

Theorem 3.2**.**

For all $z\in\mathbb{R}^{d}\setminus\overline{D}$

[TABLE]

Notice that since $|F|^{1/2}$ is a positive operator the function $W_{\text{{FDSM}}}(z)$ is strictly positive for all $z$ and as the sampling point moves away from the boundary $W_{\text{{FDSM}}}(z)$ decays. Using the eigensystem for $F$ we have that

[TABLE]

compared to the indicator given by the FM defined as

[TABLE]

Since the eigenvalues tend to zero rapidly one should avoid dividing by them in practice. The function $W_{\text{{FDSM}}}(z)$ is computed by multiplying by the square roots of the eigenvalues which will be more stable in the presence of errors in the measured far-field data. One drawback is that Theorems 2.1 and 3.2 only give that the new indicators decay as the sampling point moves away from the scatterer which may result in low contrast reconstructions. In [22] it is seen that this can be overcome by raising the new indicators to the power $p>1$ to sharpen the resolution. To recover the scatterer one can take a level curve of the indicator $W^{p}$.

Notice that the analysis in this section only requires that the far-field operator is injective with dense range and is normal which implies the orthonormal eigenvalue decomposition. The fact that the corresponding scattering operator is unitary is the key component to the analysis in this section. Therefore, this equivalence of the indicators holds for any scattering problem where the corresponding scattering operator is unitary. This is true for the inverse obstacle scattering problem where the scattered field solves

[TABLE]

along with the Sommerfeld radiation condition. Here the boundary operator is given by

[TABLE]

where $\gamma\in L^{\infty}(\partial D)$ is non-negative. This implies that $W_{\text{{FDSM}}}(z)$ can be used to recover sound soft/sound hard, isotropic, and impedance type scatterers.

4 A Tikhonov regularization based direct sampling method

The operator $|F|^{1/2}$ has been used to recover the scattering objects in previous studies where one solves the ill-posed equation

[TABLE]

which is solvable if and only if the sampling point $z\in D$. One of the main ideas proposed in [15] is to derive an approximate solution operator to the above equation and use the approximate solution operator to define a DSM. In [15] the authors approximate the solution operator using a Neumann series. Using a Neumann series to approximate the solution operator amounts to constructing a polynomial that when applied to the operator acts as the solution operator of (6). The main idea we exploit is to construct a polynomial that when evaluated at the operator $|F|$ acts as an approximate solution operator for (6). Here we propose approximating the solution operator using Tikhonov regularization, which is commonly used in the literature to solve (6) (see for e.g. [8]).

The analysis in this section again appeals to the eigenvalue decomposition. Now, recall that the orthonormal eigenvalue decomposition of the injective far-field operator $F$ is given by $(\lambda_{j},\psi_{j})\in\mathbb{C}\setminus\{0\}\times L^{2}(\mathbb{S})$. Therefore, we can define $|F|^{p}$ for $p>0$ by

[TABLE]

where the set $\{\psi_{j}\}$ is an orthonormal basis in $L^{2}(\mathbb{S})$. Note that the Tikhonov regularized solution of (6) will be denoted $g_{z}^{\alpha}$ and is the unique minimizer of the functional

[TABLE]

where $\alpha>0$ is the regularization parameter. Simple calculations give that the minimizer $g_{z}^{\alpha}$ satisfies the equation

[TABLE]

and using the eigenvalue decomposition we have that

[TABLE]

We now define the function $\Gamma_{\alpha}(t)=\frac{\sqrt{t}}{\alpha+t}$ which is continuous on the interval $[0,\|F\|]$ and we have the solution operator for the Tikhonov regularization of (6) given by the mapping

[TABLE]

In general, the regularization parameter $\alpha$ is taken to be small. The parameter $\alpha$ is chosen to be a fixed but small parameter throughout all the calculations in our experiments provided in Section 5.

In order to approximate the solution operator in (7) we exploit the fact that for all $\alpha>0$ the function $\Gamma_{\alpha}(t)$ is actually continuous for all $t\geq 0$ and therefore we have that for every $\varepsilon>0$ there is a polynomial {\color[rgb]{1,0,0}P_{\alpha,\varepsilon}(t)} such that

[TABLE]

The approximation of the solution operator is now defined by

[TABLE]

We have defined the polynomial of the operator $|F|$ via the eigenvalue decomposition as is commonly done in Linear Algebra. The Tikhonov indicator we propose in this section for a fixed $\alpha$ positive is defined as the function

[TABLE]

where {\color[rgb]{0,0,0}P_{\alpha,\varepsilon}}(t) is a polynomial. By definition of {\color[rgb]{0,0,0}P_{\alpha,\varepsilon}}(|F|)\phi_{z} we have that

[TABLE]

by appealing to the fact that $\{\psi_{j}\}$ is an orthonormal set in $L^{2}(\mathbb{S})\,.$

The goal now is to see how the new Tikhonov indicator $W_{\text{TDSM}}(z)$ compares to the indicators studied in the previous sections. To do so, we assume that the regularization parameter is known and fixed and note that for the polynomial $P_{\alpha,\varepsilon}(t)$ satisfying (8) we have that

[TABLE]

Since $\Gamma_{\alpha}(t)$ is continuous for $t\geq 0$, we let

[TABLE]

which only depends on $\alpha$ and $\|F\|$. Using basic calculus one can easily show that

[TABLE]

Now assume that $0<\varepsilon<1$ and estimate

[TABLE]

Notice that by definition $\Gamma^{2}_{\alpha}(|\lambda_{j}|)\leq|\lambda_{j}|/\alpha^{2}$ which gives that

[TABLE]

where $C(\alpha,d)$ is a positive constant depending on the regularization parameter and the dimension. By definition of $|F|^{1/2}$ it is clear that

[TABLE]

Therefore, by appealing to estimates in Section 3 we have the following result.

Theorem 4.1**.**

For all fixed $\alpha>0$ we have

[TABLE]

where the constant $\mu$ is defined in (5) and the polynomial {\color[rgb]{0,0,0}P_{\alpha,\varepsilon}}(t) satisfies (8).

The above theorem suggests that the Tikhonov indicator $W_{\text{TDSM}}(z)$ decays at least as fast as the standard DSM indicator $W_{\text{DSM}}(z)$ (up to order $\varepsilon$) and therefore we expect that $W_{\text{TDSM}}(z)$ has to approximately decay at a rate of $\mathrm{dist}(z,D)^{1-d}$ by Theorem 2.1 which implies that

[TABLE]

as $\varepsilon\to 0$. Notice that again the analysis of this section requires that the far-field operator is injective with dense range and is normal which holds for sound soft/sound hard, isotropic, and impedance type scatterers. Therefore, the above estimate for the indicator $W_{\text{TDSM}}(z)$ holds for these scatterers as well.

5 Numerical validation of the indicators

In this section, we provide some numerical examples for the three indicator functions studied in the previous sections. To this end, we will use synthetic far-field data in our numerical experiments. All of our experiments are done with MATLAB 2018a on an iMac with a 4.2 GHz Intel Core i7 processor with 8GB of memory. We will denote the discretized far-field operator as

[TABLE]

where $M$ is the number of incident and observation directions. Here the incident and observation directions are uniformly spaced on the unit circle/sphere. We give examples with random noise added to the simulated data for $u^{\infty}(\hat{x}_{i},\hat{y}_{j})$. The random noise level is given by $\delta$ where the noise is added to the matrix such that

[TABLE]

with the complex-valued random matrix $\mathbf{E}$ satisfying $\|\mathbf{E}\|_{2}=1\,.$ Here $\|\,\cdotp\|_{2}$ denotes the spectral norm of a given matrix.

To evaluate our indicators we also define the vector

[TABLE]

Therefore, we have that the DSM indicator $\left|(F\phi_{z},\phi_{z})_{L^{2}(\mathbb{S})}\right|$ is computed numerically by

[TABLE]

where the ∗ denotes the conjugate transpose and $\widetilde{W}_{\text{DSM}}(z)$ denotes the approximation of $W_{\text{DSM}}(z)$. For the indicator function based on the FM we need to numerically compute the absolute value of ${\bf F}_{\delta}$ to the half power. To do so, we let

[TABLE]

be the singular value decomposition of the matrix. Using the identity $|{\bf F}_{\delta}|^{1/2}=({\bf F}_{\delta}^{*}{\bf F}_{\delta})^{1/4}$ we have that $|{\bf F}_{\delta}|^{1/2}={\bf V}{\bf S}^{1/2}{\bf V}^{*}$. Therefore, we can numerically approximate the indicator $\left(|F|^{1/2}\phi_{z},\phi_{z}\right)_{L^{2}(\mathbb{S})}$ using the singular value decomposition such that

[TABLE]

where $(s_{j},{\bf v}_{j})\in\mathbb{R}_{>0}\times\mathbb{C}^{M}$ are the singular values and right singular vectors of ${\bf F}_{\delta}\,.$

Now for the indicator based on the Tikhonov regularization we need to construct a polynomial $P_{\alpha}(t)$ such that for all $t\in[0,\|{\bf F}_{\delta}\|]$ approximates the function $\Gamma_{\alpha}(t)$ defined in the previous section. Due to the high condition number for polynomial interpolation we notice that since $\Gamma_{\alpha}(0)=0$ we construct an approximating polynomial that has zero as a root. This gives one coefficient less to compute. In our experiments we compute {\color[rgb]{0,0,0}P_{\alpha,\varepsilon}}(t), where

[TABLE]

with $\ell=1,\ldots,10$ and $t_{\ell}$ are equally spaced point in the interval $[0,\|{\bf F}_{\delta}\|]$. One normally attempts to pick an optimal regularization parameter $\alpha$ via Morozov’s discrepancy principle. Here, we only require $\alpha$ to be positive for Theorem 4.1 to hold which is all one needs to establish the decay of the indicator function. In general, the parameter $\alpha$ is taken to be small such that $\Gamma_{\alpha}(t)\approx 1/\sqrt{t}$ holds in order to make the mapping (7) the approximate solution of (6). To this end, we fix $\alpha=10^{-2}$ ad hoc in our experiments and use a spectral cut-off to compute the coefficients $c_{k}$ where the cut-off parameter is fixed to be $10^{-8}$ in all the examples. Once {\color[rgb]{0,0,0}P_{\alpha,\varepsilon}}(t) is computed, we can numerically approximate the indicator \big{\|}{\color[rgb]{0,0,0}P_{\alpha,\varepsilon}}(|F|)\phi_{z}\big{\|}^{2}_{L^{2}(\mathbb{S})} such that

[TABLE]

where we have again used the singular values $s_{j}$ and right singular vectors ${\bf v}_{j}$. Notice that {\color[rgb]{0,0,0}P_{\alpha,\varepsilon}}(t) need only be constructed once and is continuously used to evaluate the indicator. Two interesting questions arise about the implementation of this method: how to construct the polynomial {\color[rgb]{0,0,0}P_{\alpha,\varepsilon}}(t) and how should one choose the regularization parameter? These two questions are not discussed here, but could lead to an interesting numerical and analytical investigation in the future.

Here we normalize the three indicator functions by dividing by their $L^{\infty}$ norms. In our experiments we see that the functions will be approximately one near the boundary of the scatterer and decay as the sampling point $z$ moves away from the scatterer. Therefore, one can choose a level curve to recover the boundary of the scatterer which should be taken between $0.8$ and $0.9$ for the 2D case and between $0.6$ and $0.8$ for the 3D case from our experiments. In the following subsections we see that in some cases the indicators $\widetilde{W}_{\text{FDSM}}(z)$ and $\widetilde{W}_{\text{TDSM}}(z)$ seem to give a better contrast in the reconstructions which gives a more detailed approximation of the shape of the scatterer.

5.1 Numerical results in two dimensions

Here we consider reconstructing small isotropic scatterers (i.e. $A=I$) in $\mathbb{R}^{2}$. We have that the scattered field is given by the solution to the Lippmann-Schwinger integral equation (see for e.g. [11])

[TABLE]

where $H^{(1)}_{0}$ is the first kind Hankel function of order zero and $u$ is the total field for (1)–(2). We consider the boundary of the domain $\partial D=r(\theta)\left(\cos(\theta),\sin(\theta)\right)$ where $r(\theta)$ is given by

[TABLE]

See Figure 1 for a plot of the pear-shaped, star-shaped and peanut-shaped domain.

Since the domains under consideration are small it is well known that the Born approximation, which is the first term in the Neumann series solution to the Lippmann-Schwinger integral equation, approximates the scattered field. It is given by

[TABLE]

This gives that the far-field pattern can be approximated by

[TABLE]

which is a good approximation for the far-field pattern since $|D|\ll 1$. To evaluate the integral we use the built-in numerical 2D integrator ‘integral2’ in MATLAB.

In Figures 2–4 we fix the wave number $k=10$ and the refractive index $n=1/2$. We add $5\%$ random noise to the approximated far-field pattern which corresponds to $\delta=0.05$. There are $M=32$ uniformly spaced incident and observation directions given by $\hat{y}_{j}=\hat{x}_{j}=\left(\cos(\theta_{j}),\sin(\theta_{j})\right)$ where $\theta_{j}$ are uniformly spaced points in $[0,2\pi)$. The sampling region is given by $[-1,1]\times[-1,1]$ where the sampling points are taken to be $100\times 100$ equally spaced points in the sampling region. Computing the three indicator functions takes roughly three seconds for each of the domains under consideration. In Figures 2–4 we plot the three indicator functions studied in the previous section. We see in Figures 2–4 that the two indicators $\widetilde{W}_{\text{FDSM}}(z)$ and $\widetilde{W}_{\text{TDSM}}(z)$ give comparable results to what is given by $\widetilde{W}_{\text{DSM}}(z)$.

We additionally show the robustness of the newly proposed indicator based on the Tikhonov regularization for the peanut-shaped scatterer. We use the noise levels $\delta=1\%$, $5\%$, and $10\%$ for the far-field data and show the results for the indicator approximation $\widetilde{W}_{\text{TDSM}}(z)$ in Figure 5. As we can see, the reconstructions of the peanut-shaped scatterer are robust with respect to these noise levels. Almost no difference can be seen visually. It is well-known that the robustness is tied to the choice of the regularization parameter which is not studied here. However, in practice one would like to pick an optimal $\alpha$ depending on the noise level $\delta$ which is left for future investigation.

5.2 Numerical results in three dimensions

We now give some numerical examples in $\mathbb{R}^{3}$ for the case of $A=\alpha I$ where both $\alpha$ and $n$ are constants. In order to compute the far-field pattern we use boundary integral equations derived from Green’s representation formula for the scattered field $u^{s}$ in $\mathbb{R}^{3}\setminus\overline{D}$ and the total field $u$ in $D$ (see [1] for details). This gives a 2$\times$2 system of boundary integral equations on the boundary $\partial D$ for the densities $\varphi_{1}$ and $\varphi_{2}$. We employ a boundary element collocation method to solve the system numerically. Here, the far-field operator is approximated for simplicity by constant interpolation over a triangulation of the unit sphere. For a possible higher-order approximation of it, we refer the reader to [20].

In our numerical examples the computed densities $\varphi_{1}$ and $\varphi_{2}$ depend on the incident direction $\hat{y}$ since each incident direction corresponds to a new right-hand side for the boundary integral equations. The corresponding far-field pattern $u^{\infty}$ can be derived from Green’s representation and is given by

[TABLE]

The boundary of the domain is given in spherical coordinates such that

[TABLE]

where $r(\phi)$ is given by $r(\phi)=1$,

[TABLE]

which represent a sphere, peanut-shaped, and acorn-shaped domain, respectively. See Figure 6 for a triangularization of their surfaces. Note that some of the far-field data have already been used in [3, Section 8].

In Figures 7–9 we fix the wave number $k=2$ with $n=1/2$ and $A=2I$ in our examples. We add $5\%$ random noise to the approximated far-field pattern which corresponds to $\delta=0.05$. Here, we use $M=258$ incidence and observation directions that are ‘almost’ uniformly spaced on the unit sphere (see [17, Appendix A.1] for details). We plot the three indicators on the $y-z$ plane such that the sampling region is given by $[-2,2]\times[-2,2]$ where the sampling points are taken to be $100\times 100$ equally spaced points in the sampling region. Computing the three indicator functions takes roughly five to six seconds for each of the domains under consideration. In Figures 7–9 we contour plots of the three indicator functions for spherical, peanut-shaped, and acorn-shaped domain. Again we see that the three indicators give comparable reconstructions for the sphere and peanut-shaped obstacle. In Figure 9 we see that the indicator $\widetilde{W}_{\text{TDSM}}(z)$ gives a more detailed reconstruction of the acorn-shaped scatterer.

Finally, we show the reconstructions for the peanut-shaped domain by the newly proposed indicator based on the Tikhonov regularization using the noise levels $\delta=1\%$, $5\%$, and $10\%$ for the far-field data. The results for the indicator approximation $\widetilde{W}_{\text{TDSM}}(z)$ are given in Figure 10. Here, we can see that we obtain the best reconstruction for the noise level $1\%$. The reconstruction deteriorates significantly as the noise level increases.

6 Summary and outlook

In this article, we have studied three DSM indicator functions using the theoretical basis of the FM for anisotropic materials with real-valued coefficients. The equivalence of two indicators previously studied is shown as well as a new indicator based on the Tikhonov regularization applied to the far-field equation for the factorization method is proposed. Precisely, we are able to prove that one of the indicator functions decays as the sampling point moves away from the scatterer. Note that the results here are stated for anisotropic materials but the analysis is valid whenever the scattering operator is unitary and the far-field operator is injective with dense range which means that these results hold for a wide range of scattering objects. Numerical examples show the validity of the three indicator functions. In our experiments we are able to reconstruct the scatterers in seconds which gives that these DSMs are computationally cheap and rigorously justified. The analysis in this article depends mainly on the orthonormal spectral decomposition of the far-field operator. It is known that for the case of complex-valued coefficients the far-field operator fails to be normal and therefore does not have the orthonormal spectral decomposition which is vital in our approach. However, one can still apply the FM either to the operators $\mathrm{Im}(F)$ or $F_{\sharp}=|\mathrm{Re}(F)|+|\mathrm{Im}(F)|$ (see for e.g. [8, 17, 18, 21]). By definition one has that these operators are self-adjoint and compact which implies they have an orthonormal spectral decomposition. One can then construct DSM indicator functions using either $\mathrm{Im}(F)$ or $F_{\sharp}$ also for the case of scatterers with complex-valued coefficients.