Imaging point sources in heterogeneous environments

TL;DR

This paper studies how uncertainties in the environment's refractive index affect the accuracy of imaging point sources using boundary or far-field measurements, providing stability results and numerical validation.

Contribution

It extends existing imaging theories by analyzing the impact of environment uncertainty on source reconstruction stability in heterogeneous media.

Findings

Reconstruction stability is maintained under measurement errors.

Smooth variations in the environment have limited impact on imaging accuracy.

Numerical simulations confirm theoretical stability results.

Abstract

Imaging point sources in heterogeneous environments from boundary or far-field measurements has been extensively studied in the past. In most existing results, the environment, represented by the refractive index function in the model equation, is assumed known in the imaging process. In this work, we investigate the impact of environment uncertainty on the reconstruction of point sources inside it. Following the techniques developed by El Badia and El Hajj (C. R. Acad. Sci. Paris, Ser. I, 350 (2012), 1031-1035), we derive stability of reconstructing point sources in heterogeneous media with respect to measurement error as well as smooth changes in the environment, that is, the refractive index. Numerical simulations with synthetic data are presented to further explore the derived stability properties.

Figures

Tables

Table 1. Table 1: True and reconstructed intensities of the point sources, ( λ 1 , λ 2 , λ 3 , λ 4 ) subscript 𝜆 1 subscript 𝜆 2 subscript 𝜆 3 subscript 𝜆 4 (\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}) , in different numerical experiments.

\@tabular@row@before@xcolor  \@xcolor@tabular@before Experiment	True value	Reconstructions with noisy data
\@tabular@row@before@xcolor  \@xcolor@row@after		1% noise	5% noise
\@tabular@row@before@xcolor  \@xcolor@row@after1: medium (33)	(0.89,0.73,0.71,0.52)	(0.91,0.73,0.68,0.54)	(0.92,0.65,0.76,0.46)
\@tabular@row@before@xcolor  \@xcolor@row@after1: medium (34)	(0.89,0.73,0.71,0.52)	(0.88,0.72,0.71,0.52)	(0.88,0.72,0.69,0.54)
\@tabular@row@before@xcolor  \@xcolor@row@after2	(0.89,0.73,0.71,0.52)	(0.91,0.76,0.65,0.58)	(0.87,0.83,0.56,0.64)
\@tabular@row@before@xcolor  \@xcolor@row@after3	(0.89,0.73,0.71,0.52)	(0.90,0.81,0.65,0.56)	(0.89,0.85,0.61,0.60)
\@tabular@row@before@xcolor  \@xcolor@row@after
\@tabular@row@before@xcolor  \@xcolor@row@after The simulations show that one can indeed reconstruct point sources, both their locations and their intensities, inside heterogeneous media when the media are not unreasonably complex. We want to emphasize here that in our theoretical analysis as well as numerical simulations, both the true sources and the sources to be reconstructed are explicitly assumed to be point sources. In other words, we explicitly search for the locations and strengths of the point sources, instead of reconstructing spatially distributed sources hoping that the result will give us point sources. In general, we observe from our extensive numerical simulations that when the refractive index is exactly known, the reconstructions are quite stable when the number of point sources is small. However, the reconstructions become too sensitive to algorithmic parameters when the number of point sources gets large. We also want to emphasize that it is important to impose the constraints on the separability of the point sources in the numerical simulations. In other words, we have to explicitly ensure that the point sources to be reconstructed are far away from each other. Even in this case, the reconstructions are sensitive to the initial guess of the locations of the point sources. The objective function that we minimize to reconstruct the point sources can not differentiate between the true point sources and the equivalent class of re-labeled point sources. Therefore, the minimization algorithm could be easily fooled to jump between different intermediate configurations if the point sources are not well-separated. To further illustrate on this issue, we plot in Figure 3 the (normalized) objective functional $\Psi$ defined in (28) as a function of the location of a single point source (the intensity of the source being assumed known). The true point source is located at $(0.443,0.298)$. While it is clear from the plot that the objective function is convex with respect to the location of the point source when $k=5$, this is not true anymore when $k=8$ and $k=12$. In the case of $k=8$, two local minmizers emerge at $y\approx 0.85$. More local minimizers emerge when $k=12$. These plots show that even in the case of a single point source, when the initial guess is far from the true position, the minimization algorithm could return wrong reconstructions. We can not visualize this phenomenon in the case of more than one point source. However, one can easily imagine that the situation would be far worse in that scenario. Figure 3: Normalized objective function $\Psi$, defined in  (28), as a function of the location of a single point source. The true location is at $(0.443,0.298)$. Shown from left to right are the three cases of $k=5$, $k=8$ and $k=12$ respectively. Experiment 2 [Recovery in an Unknown Environment]. In the second set of simulations, we reconstruct point sources in medium (33) assuming that both the medium and the point sources are not known. The reconstructions of the medium and the locations of the sources are shown in Figure 3 and the reconstructed intensities of the point sources are summarized in the third row of Table 3.
\@tabular@row@before@xcolor  \@xcolor@row@after Let us emphasize again that the reconstruction here is done in two steps. In the first step, we reconstruct the refractive index using multiple differential data sets. In the second step, we fix the refractive index, which is the reconstructed one, and reconstruct the point sources from one Cauchy data set. If we compare the reconstructions in Figure 3 with those in Figure 1 (which are reconstructed under the true medium $n$ in (33)), and the reconstructions of intensities in the third row of Table. 3 with those in the first row of the same table, we see that the reconstructions of the point sources are different but are of similar quality. That is, smooth changes in the refractive index introduces relatively small error in the reconstruction of point sources. This confirms our stability result in Theorem 2.6. Experiment 3 [Recovery in an Unknown Environment]. We repeat here the simulations in Experiment 2 for the medium (34). The reconstructions of the medium and the locations of the sources are shown in Figure LABEL:FIG:Unknown_Media_2 and the reconstructed intensities of the point sources are summarized in the last row of Table 3.

Equations

$$\begin{array}[]{rcll}\Delta u+k^{2}\big{(}1+n(\mathbf{x})\big{)}u&=&q(\mathbf{x}),&\mbox{in}\ \ \Omega\\ u&=&f(\mathbf{x}),&\mbox{on}\ \partial\Omega\end{array}$$

$$q(x)=\sum _{j=1}^m{\lambda }_j\delta (\mathbf{x}-{\mathbf{x}}_j).$$

$$g(\mathbf{x})=\left({\dfrac{\partial u}{\partial\nu}}\right){\bigg{|}}_{\partial\Omega}\equiv{\boldsymbol{\nu}}\cdot\nabla u|_{\partial\Omega},$$

$$q_{\ell }=\sum _{j=1}^{m_{\ell }}{\lambda }_{\ell ,j}\delta (\mathbf{x}-{\mathbf{x}}_{\ell ,j}),\ell =1,2$$

$${\mathbf{x}}_{1,j}={\mathbf{x}}_{2,\pi (j)},{\lambda }_{1,j}={\lambda }_{2,\pi (j)},1\le j\le m$$

$$\begin{array}[]{cl}\Delta w+k^{2}(1+n)w=0,&\mbox{in}\ \Omega_{\varepsilon}:=\Omega\backslash\bigcup_{\ell,j}D_{\varepsilon}^{\ell,j}\\ \\ w=0,\quad\partial_{\nu}w=0,&\mbox{on}\ \partial\Omega\end{array}$$

$$\Delta \varphi +k^2(1+n)\varphi =0,\text{mbox}in\Omega$$

$$\underset{\mathbf{x}\in \Omega }{\sup }(|\varphi (\mathbf{x})|+|\nabla \varphi (\mathbf{x})|)<\mu \underset{\mathbf{x}\in \Omega }{\inf }|\varphi (\mathbf{x})|.$$

$$\varphi (\mathbf{x})=e^{i\zeta \cdot \mathbf{x}}(1+r(\mathbf{x})),\text{ with }\Vert r{\Vert }_{{\mathcal{H}}^3(\Omega )}\le \mathfrak{c}\frac{k^2}{|\zeta |}\Vert n{\Vert }_{{\mathcal{H}}^3(\Omega )}$$

$$\underset{\mathbf{x}\in \Omega }{\sup }(|\varphi (\mathbf{x})|+|\nabla \varphi (\mathbf{x})|)<\mu \underset{\mathbf{x}\in \Omega }{\inf }|\varphi (\mathbf{x})|.$$

$$\mathcal{S}(\mathbf{x})=\overline{({\mathfrak{u}}_2+\mathfrak{i}{\mathfrak{u}}_3)}.$$

$$\Delta \psi +k^2(1+n)\psi =0,\text{mbox}in\Omega$$

$$\psi ({\mathbf{x}}_j)=0,1\le j\le m.$$

$$\psi(\mathbf{x})=\phi(\mathbf{x})\prod_{j=1}^{m}\big{(}\mathcal{S}(\mathbf{x})-\mathcal{S}(\mathbf{x}_{j})\big{)}.$$

$$\begin{gathered} {\rho }_{\mathbf{x}}:=\underset{1\le j\le m}{\max }|\mathcal{S}({\mathbf{x}}_{1,j})-\mathcal{S}({\mathbf{x}}_{2,\pi (j)})| \\ \le {\mathfrak{c}}_1{\left(\frac{\sqrt{|\partial \Omega |}({\mathrm{diam}}_{\mathcal{S}}(\Omega ){)}^{2m-1}}{\underline{\lambda }{\sigma }^{m-1}}\left(\Vert g_1-g_2{\Vert }_{L^2(\partial \Omega )}+\Vert f_1-f_2{\Vert }_{L^2(\partial \Omega )}\right)\right)}^{\frac{1}{m}} \end{gathered}$$

$$\rho_{\lambda}:=\max_{1\leq j\leq m}|\lambda_{1,j}-\lambda_{2,\pi(j)}|\\ \leq\mathfrak{c}_{2}\overline{\lambda}\widetilde{\rho}_{\mathbf{x}}+\mathfrak{c}_{3}\sqrt{|\partial\Omega|}({\rm diam}_{\mathcal{S}}(\Omega))^{2m-2}\big{(}\|g_{1}-g_{2}\|_{L^{2}(\partial\Omega)}+\|f_{1}-f_{2}\|_{L^{2}(\partial\Omega)}\big{)}.$$

$${\psi }_{j^{\prime }}(\mathbf{x})=\varphi (\mathbf{x})\prod _{i=1}^m(\mathcal{S}(\mathbf{x})-\mathcal{S}({\mathbf{x}}_{1,i}))\prod _{j=1,j\ne j^{\prime }}^m(\mathcal{S}(\mathbf{x})-\mathcal{S}({\mathbf{x}}_{2,j}))$$

$$\Delta {\psi }_{j^{\prime }}+k^2(1+n){\psi }_{j^{\prime }}=0,\text{mbox}in\Omega$$

$${\psi }_{j^{\prime }}({\mathbf{x}}_{\ell ,j})=0,\forall (\ell ,j)\ne (2,j^{\prime }).$$

$$\displaystyle\left|\int_{\partial\Omega}\big{(}(f_{1}-f_{2})\partial_{\nu}\psi_{j^{\prime}}-\psi_{j^{\prime}}(g_{1}-g_{2})\big{)}d\mathfrak{s}\right|=\left|\lambda_{2,j^{\prime}}\psi_{j^{\prime}}(\mathbf{x}_{2,j^{\prime}})\right|$$

$$=$$

$$\left|\int_{\partial\Omega}\big{(}(f_{1}-f_{2})\partial_{\nu}\psi_{j^{\prime}}-\psi_{j^{\prime}}(g_{1}-g_{2})\big{)}d\mathfrak{s}\right|\\ \leq\|g_{1}-g_{2}\|_{L^{2}(\partial\Omega)}\|\psi_{j^{\prime}}\|_{L^{2}(\partial\Omega)}+\|f_{1}-f_{2}\|_{L^{2}(\partial\Omega)}\|\partial_{\nu}\psi_{j^{\prime}}\|_{L^{2}(\partial\Omega)},$$

$$\Vert {\psi }_{j^{\prime }}{\Vert }_{L^2(\partial \Omega )}$$

$$\Vert {\partial }_{\nu }{\psi }_{j^{\prime }}{\Vert }_{L^2(\partial \Omega )}$$

$$\underline{\lambda}\rho_{j^{\prime}}^{m}\sigma^{m-1}\leq\widetilde{\mathfrak{c}}\sqrt{|\partial\Omega|}({\rm diam}_{\mathcal{S}}(\Omega))^{2m-1}\big{(}\|g_{1}-g_{2}\|_{L^{2}(\partial\Omega)}+\|f_{1}-f_{2}\|_{L^{2}(\partial\Omega)}\big{)}$$

$$\underset{1\le j^{\prime }\le m}{\max }{\rho }_{j^{\prime }}\le {\left(\frac{\tilde{\mathfrak{c}}\sqrt{|\partial \Omega |}({\mathrm{diam}}_{\mathcal{S}}(\Omega ){)}^{2m-1}}{\underline{\lambda }{\sigma }^{m-1}}\left(\Vert g_1-g_2{\Vert }_{L^2(\partial \Omega )}+\Vert f_1-f_2{\Vert }_{L^2(\partial \Omega )}\right)\right)}^{1/m}.$$

$$\max \{\underset{1\le j^{\prime }\le m}{\max }\underset{1\le j\le m}{\min }|\mathcal{S}({\mathbf{x}}_{2,j^{\prime }})-\mathcal{S}({\mathbf{x}}_{1,j})|,\underset{1\le j^{\prime }\le m}{\max }\underset{1\le j\le m}{\min }|\mathcal{S}({\mathbf{x}}_{1,j^{\prime }})-\mathcal{S}({\mathbf{x}}_{2,j})|\}.$$

$${\phi }_{j^{\prime }}(\mathbf{x})=\varphi (\mathbf{x})\prod _{j=1,j\ne j^{\prime }}^m(\mathcal{S}(\mathbf{x})-\mathcal{S}({\mathbf{x}}_{1,j}))(\mathcal{S}(\mathbf{x})-\mathcal{S}({\mathbf{x}}_{2,\pi (j)})).$$

$$\Delta {\phi }_{j^{\prime }}+k^2(1+n){\phi }_{j^{\prime }}=0,\text{mbox}in\Omega$$

$${\phi }_{j^{\prime }}({\mathbf{x}}_{\ell ,j})=0,\forall (\ell ,j)\notin \{(1,j^{\prime }),(2,\pi (j^{\prime }))\}.$$

$$\left|{\lambda }_{2,\pi (j^{\prime })}{\phi }_{j^{\prime }}({\mathbf{x}}_{2,\pi (j^{\prime })})-{\lambda }_{1,j^{\prime }}{\phi }_{j^{\prime }}({\mathbf{x}}_{1,j^{\prime }})\right|=\left|{\int }_{\partial \Omega }\left((f_1-f_2){\partial }_{\nu }{\phi }_{j^{\prime }}-{\phi }_{j^{\prime }}(g_1-g_2)\right)d\mathfrak{s}\right|.$$

Full text

Key words. Inverse source problems, Helmholtz equation, point sources, stability estimates, numerical reconstructions, uncertainty characterization, inverse problems

AMS subject classifications 2010. 15A29, 35R30, 49N45, 65N21, 78A46.

1 Introduction

Recovering radiative sources inside heterogeneous media from boundary or far-field measurements has applications in many branches of science and technology [AmBaFl-SIAM02, Baltes-Book78, BaLiTr-JDE10, BaBeBeLe-IP05, CaEw-IPBVP75, Chin-ACM10, ElVa-IP09, GaOsTa-JIIPP13, HoCl-MB97, Ikehata-IP99, ImYa-IP98, IsLu-IPI18, JiLiYa-JDE17, Larsen-JQSRT73, MaTs-IP13, Nicaise-SIAM00, PuYa-IP96, Sezginer-IP87, Siewart-JQSRT93, Yamamoto-IP95]. Extensive mathematical and computational studies of such inverse source problems have been performed in the past decades; see, for instance,  [AnBuEr-Book97, Isakov-Book90, Isakov-Book06] and references therein for recent reviews on the subject.

In this work, we are interested in a source recovery problem where the source to be reconstructed is the superposition of point sources [BaBeBeLe-IP05, CaEw-IPBVP75, ElBadia-IP05, ElHa-JIIPP02, ElEl-CRASP12, ElNa-IP11, FaHaEs-IPI13, GaZhCoWa-BOE10, KaLe-IP04, KoYa-IP02, LiTa-CiCP09, MaTs-IP13, OhInOh-IP11, Vessella-IP92]. Unlike general source functions, point sources are efficiently characterized by their locations and strengths, a fact that significantly reduces the dimension of the parameter space of the inverse problems. This dimension reduction often enables one to obtain uniqueness in the inverse problem with minimum amount of observed data and provides the possibility of utilizing efficient reconstruction algorithms, for instance these based on compressive sensing [ChMoPa-IP12, FaStYa-SIAM10], in the source recovery process.

To formulate our problem, let $\Omega\subset\mathbb{R}^{d}$ ($d\geq 3$) be a simply connected domain with $\mathcal{C}^{2}$ boundary $\partial\Omega$. Let $u(\mathbf{x})$ be the solution to the following boundary value problem to the Helmholtz equation:

[TABLE]

where the real-valued function $n(\mathbf{x})$ is the refractive index, $q(\mathbf{x})$ and $f(\mathbf{x})$ are internal and boundary source terms respectively. We assume that $n(\mathbf{x})$ has a compact support in $\Omega$, that is, $\operatorname{supp}(n)\subset\subset\Omega$, and that $1+n(\mathbf{x})>0$, $\forall\mathbf{x}\in\Omega$. We assume that [math] is not an eigenvalue of the operator $\Delta+k^{2}(1+n)$ with homogeneous Dirichlet boundary condition such that the problem (1) admits a unique solution for given source functions $q$ and $f$.

We assume that the internal source function $q(\mathbf{x})$ is a superposition of $m$ point sources located at $\{\mathbf{x}_{j}\}_{j=1}^{m}$ with strengths $\{\lambda_{j}\}_{j=1}^{m}$, that is,

[TABLE]

The strengths $\{\lambda_{j}\}$ are all assumed to be real-valued so that there is no physical absorption occurring at the point sources.

The Helmholtz equation (1) can be viewed as a simplified frequency-domain model for either electromagnetic or ultrasound wave propagation, depending on the value of the parameters, mainly the wavenumber $k$, in the equation. The mathematical derivations in the rest of the paper implicitly assume that the wavenumber $k$ is real-valued and $k>0$. We believe, however, the same types of calculations can be carried out for the zero-frequency case ($k=0$) as those studied in  [BaBeBeLe-IP05, ElBadia-IP05] or in [WaLiJi-MP05] with an extra absorption term.

Let us also mention that since linearizing inverse coefficient problems often results in inverse source problems, the point source reconstruction problem we study in this paper is closely related to the problem of reconstructing small volume inclusions in background media [AmKa-IP03, AmMoVo-ESAIM03] and the problem of imaging small scatterers in complex media [BoPaTs-IP03, BoPaTs-IP05, ChMoPa-IP16, DeMaGr-JASA05]. The main difference is that in our derivation below, we can utilize explicitly the fact that point sources are singular. The secondary sources created by small scatterers or inclusions, however, do not carry the same level of singularity of the point sources.

We are interested in the problem of reconstructing the point sources, i.e. their locations and strengths, from Cauchy data $(f,g)$ where the boundary measurement $g(\mathbf{x})$ is given as

[TABLE]

${\boldsymbol{\nu}}(\mathbf{x})$ being the unit outer normal vector of the domain boundary at $\mathbf{x}\in\partial\Omega$.

To our best knowledge, in all the previous work on point source recovery, the environment, that is, the refractive index $n(\mathbf{x})$ in our formulation, in which point sources (or point source like localized objects), are to be sought is assumed to be known exactly, with the only exception in  [BuLaTaTs-AFR09] where the authors tried to reconstruct point sources in homogeneous media with unknown impenetrable obstacles. In other words, in the environment under which the data are collected is the same as the environment that data are back-propagated to reconstruct the point sources. Moreover, the mathematical works where stability of source reconstructions are derived with respect to noise in measured data are all done under the assumption that the environment is homogeneous (so that one could have access to the explicit form of the associated Green’s function).

In the rest of the paper, we prove a stability result, following the techniques developed by El Badia and El Hajj (C. R. Acad. Sci. Paris, Ser. I, 350 (2012), 1031-1035) [ElEl-CRASP12], on the recovering of point sources in smooth inhomogeneous environment, that is, when the refractive index $n(\mathbf{x})$ varies smoothly in space. We also prove the stability of point source reconstruction with respect to smooth changes in the environment: if we perform reconstructions in a medium that is only slightly, in appropriate sense, different than the medium from which we collected the data, then the reconstructions are only slightly different from the reconstructions in the exact medium that generated the data.

2 Main results

We now present the main result of this short paper: (i) the stability of reconstructing point sources in an heterogeneous medium; and (ii) the stability of the reconstructions with respect to smooth changes of the medium.

2.1 Stability in heterogeneous media

We first consider the case where the underlying medium is known but heterogeneous, that is, we have a known but spatially varying refractive index $n(\mathbf{x})$. Under this circumstance, we can uniquely reconstruct the source $q(\mathbf{x})$ from a single pair of Cauchy data $(f,g)$ with a given Hölder type of stability. This is a generalization of the stability results of El Badia and El Hajj established in [ElEl-CRASP12, ElNa-IP11].

We make the following general assumptions on the setup of the problem.

(A). The domain $\Omega$ has a $\mathcal{C}^{2}$ boundary $\partial\Omega$. The refractive index $n(\mathbf{x})$ is real-valued and smooth, with support $\operatorname{supp}(n)\subset\Omega$. The point sources are well separated in the sense that $\displaystyle\min_{i\neq j}|\mathbf{x}_{i}-\mathbf{x}_{j}|\geq\mathfrak{c}>0$ for some $\mathfrak{c}$. The point sources are sufficiently far away from the boundary of the domain such that ${\rm dist}(\mathbf{x}_{j},\partial\Omega)\geq\widetilde{\mathfrak{c}}>0$, $\forall j$, for some $\widetilde{\mathfrak{c}}$. The strengths of the point sources satisfy $0<\underline{\lambda}\leq\lambda_{j}\leq\overline{\lambda}<+\infty$, $\forall j$, for some $\underline{\lambda}$ and $\overline{\lambda}$. The illumination boundary source $f$ is the restriction of a $\mathcal{C}^{\infty}$ function to $\partial\Omega$.

It will be clear that the smoothness assumptions on the refractive index $n(\mathbf{x})$ and the boundary source $f$ are not completely necessary. In fact, being $\mathcal{C}^{3}$ is sufficient for all the results to hold. It should also be noted that the assumptions on the point sources imply that $\mathbf{x}_{i}\neq\mathbf{x}_{j}$ whenever $i\neq j$, a fact that is implicitly used later when we study uniqueness of reconstructions.

Lemma 2.1**.**

Let

[TABLE]

be two sets of point sources satisfying the assumptions in (A), $u_{\ell}$ ($\ell=1,2$) the corresponding solutions to the Helmholtz equation (1). Then $(f_{1},g_{1})=(f_{2},g_{2})$ implies that $m_{1}=m_{2}$ $(:=m)$ and

[TABLE]

for some permutation $\pi\in{\rm Sym(m)}$.

Proof.

This result follows from the unique continuation principle for Cauchy problems of elliptic equations. Let $D_{\varepsilon}^{\ell,j}:=\mathcal{B}(\mathbf{x}_{\ell,j},\varepsilon)$ be the disk centered at $\mathbf{x}_{\ell,j}$ with radius $\varepsilon>0$ small enough such that $D_{\varepsilon}^{\ell,j}\subset\Omega$. We define $w=u_{1}-u_{2}$ and verify that $w$ solves

[TABLE]

Since $\partial\Omega\subset\partial\Omega_{\varepsilon}$, we then conclude from the unique continuation principle [Isakov-Book06, Theorem 3.3.1] that $w(\mathbf{x})=0$, $\forall\mathbf{x}\in\Omega_{\varepsilon}$ with any $\varepsilon>0$. If we take $\varepsilon\to 0$, this implies that $w=0$ except at the locations $\mathbf{x}_{\ell,j}$. Therefore $w$ must be a finite linear combination of point sources and their derivatives. This is impossible. Therefore $q_{1}=q_{2}$ up to a possible permutation $\pi$, that is, renumbering of the point sources. ∎

We now study the stability of the reconstruction. Following [ElEl-CRASP12, ElNa-IP11], we look at the stability issue for an algebraic reconstruction technique that is based on the projection of the point sources into planes in $\mathbb{R}^{3}$. Due to the fact that the medium is heterogeneous, we need to find good ways to do the projection. In the next two lemmas, we introduce our method of projection onto surfaces determined by the solutions of the Helmholtz equation (which are controlled by the medium).

Lemma 2.2**.**

Under the assumptions in (A) on the refractive index and $\Omega$, there exists a complex-valued function $\phi(\mathbf{x})$ and a constant $\mu>0$ such that $\phi(\mathbf{x})$ solves

[TABLE]

$|\nabla\phi(\mathbf{x})|\neq 0$, $\forall\mathbf{x}\in\Omega$, and

[TABLE]

Proof.

With the regularity of $n(\mathbf{x})$ assumed in (A), we can take $\phi$ as the well-known complex geometrical optics (CGO) solution to (6); see for instance [NaUhWa-IP13]. More precisely, let $\zeta=\eta+\mathfrak{i}\xi$ with $\eta\in\mathbb{R}^{3}$ and $\xi\in\mathbb{R}^{3}$ given vectors such that $\eta\cdot\xi=0$ and $\zeta\cdot\zeta=k^{2}$ (i.e. $|\eta|^{2}-|\xi|^{2}=k^{2}$). It is shown in [NaUhWa-IP13] that (6) has a solution of the form

[TABLE]

when $|\zeta|$ is sufficiently large. Then by the Sobolev embedding theorem [AdFo-Book03], $r\in\mathcal{C}^{1,1/2}(\Omega)$. If we choose $|\zeta|$ large such that $\|r\|_{\mathcal{C}^{1,1/2}(\Omega)}<1/2$, then we will have $|\nabla\phi|\neq 0$ for all $\mathbf{x}\in\Omega$ and $\phi\in\mathcal{C}^{1,1/2}(\Omega)$. Therefore we can find a constant $\mu>0$ such that

[TABLE]

This completes the proof. ∎

Definition 2.3**.**

The function $\phi$ introduced in Lemma 2.2 defines a local frame $(\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3})$ on $\Omega$, where $\mathbf{e}_{1}=\frac{\nabla\phi}{|\nabla\phi|}\in\mathbb{C}^{3}$, and $[\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3}]$ forms an unitary matrix, which also determines a local coordinate change $\mathcal{U}:\mathbb{R}^{3}\rightarrow\mathbb{C}^{3}$, denoted by $\mathcal{U}(\mathbf{x})=(\mathfrak{u}_{1}(\mathbf{x}),\mathfrak{u}_{2}(\mathbf{x}),\mathfrak{u}_{3}(\mathbf{x}))$, from the Cartesian coordinate to the local frame. We define a projection $\mathcal{S}:\mathbb{C}^{3}\rightarrow\mathbb{C}$ as:

[TABLE]

This projection defines a pseudo distance function $\text{\rm dist}_{\mathcal{S}}(\mathbf{x},\mathbf{y})=|\mathcal{S}(\mathbf{x})-\mathcal{S}(\mathbf{y})|$ and the diameter of $\Omega$ under the projection $\mathcal{S}$ is denoted by ${\rm diam}_{\mathcal{S}}(\Omega):=\displaystyle\sup_{\mathbf{x},\mathbf{y}\in\Omega}{\rm dist}_{\mathcal{S}}(\mathbf{x},\mathbf{y})$.

The project $\mathcal{S}$ we defined here is clearly not unique in the sense that one can rotate the coordinates to use $(\mathfrak{u}_{3},\mathfrak{u}_{1})$ or $(\mathfrak{u}_{1},\mathfrak{u}_{2})$ to replace $\mathfrak{u}_{2}$ and $\mathfrak{u}_{3}$. Moreover, since $\phi$ can be chosen differently using the complex vector $\zeta$ (which controls the boundary condition needed), we can construct the $\mathcal{S}$ that we need by selection a specific vector.

The projection we just introduced allows us to construct the following test functions.

Lemma 2.4**.**

Let $\{\mathbf{x}_{j}\}_{j=1}^{m}\subset\Omega$ be arbitrary distinct points. Then under assumption (A), there exists a function $\psi(\mathbf{x})$ solving

[TABLE]

such that

[TABLE]

Proof.

We construct the function $\psi$ as follows:

[TABLE]

Then clearly $\psi(\mathbf{x}_{j})=0$, $j=1,2,\cdots,m$. It is straightforward to verify that $\Delta\mathcal{S}=0$ and $\nabla\phi\cdot\nabla\mathcal{S}=0$, which allow us to check that $\psi$ solves  (8). ∎

We are now ready to prove the stability of the reconstruction. Our reconstruction scheme follows a two-step process. In the first step, the locations of the point sources are probed by a projection method. In the second step, we use the reconstructed locations $\{\mathbf{x}_{j}\}_{j=1}^{m}$ to reconstruct the strengths of the point sources $\{\lambda_{j}\}_{j=1}^{m}$.

Theorem 2.5**.**

Let $q_{1}$ and $q_{2}$ be two sources of the form (4) with $m:=m_{1}=m_{2}$ that are reconstructed from the Cauchy data $(f_{1},g_{1})$ and $(f_{2},g_{2})$ respectively. Let $\mathcal{S}$ be the projection in Definition 2.3. Let $\sigma:=\displaystyle\min_{\ell,i\neq j}|\mathcal{S}(\mathbf{x}_{\ell,i})-\mathcal{S}(\mathbf{x}_{\ell,j})|$ and assume that $\sigma>0$. Then, under the assumptions in (A), there exists a permutation $\pi\in{\rm Sym(m)}$ acting on $\{1,2,\cdots,m\}$ and a constant $\mathfrak{c}_{1}$ depending on $\Omega$, $\phi$, $\mathcal{S}$, and $m$, such that

[TABLE]

where $|\partial\Omega|$ is the surface measure of $\partial\Omega$. Assume further that $\rho_{\mathbf{x}}<\sigma$, and let $\widetilde{\rho}_{x}:=\displaystyle\max_{1\leq j\leq m}|\mathcal{U}(x_{1,j})-\mathcal{U}(\mathbf{x}_{2,\pi(j)})|$, then there exists constants $\mathfrak{c}_{2}$ and $\mathfrak{c}_{3}$, again depending on $\Omega$, $\phi$, $\mathcal{S}$, and $m$, such that

[TABLE]

Proof.

Let $u_{\ell}$ $(\ell=1,2)$ be the solution to the Helmholtz equation (1) with source $q_{\ell}$. We define $w:=u_{1}-u_{2}$. Then $w$ solves the Helmholtz equation (1) with boundary data $(w,\partial_{\nu}w):=(f_{1}-f_{2},g_{1}-g_{2})$. Let $1\leq j^{\prime}\leq m$ be an integer. From Lemma 2.4, we find a function, $\phi$ being defined in Lemma 2.2,

[TABLE]

that solves the equation

[TABLE]

and satisfies

[TABLE]

Multiplying the equation for $w$ by $\psi_{j^{\prime}}$ and the equation for $\psi_{j^{\prime}}$ by $w$, taking the difference of the results, and applying Green’s identity, we have, with $\mathfrak{s}(\mathbf{x})$ the Lebesgue measure on $\partial\Omega$,

[TABLE]

where $\theta:=\inf_{\mathbf{x}\in\Omega}|\phi(\mathbf{x})|$, $\rho_{j^{\prime}}:=\displaystyle\min_{1\leq j\leq m}|\mathcal{S}(\mathbf{x}_{2,j^{\prime}})-\mathcal{S}(\mathbf{x}_{1,j})|$, and we have used $|\lambda_{2,j^{\prime}}|\geq\underline{\lambda}$ by the assumptions in (A).

Meanwhile, by the Cauchy-Schwartz inequality, we have

[TABLE]

where we can estimate

[TABLE]

with $\vartheta:=\|\phi\|_{W^{1,\infty}(\Omega)}$ and $\mathfrak{c}=\mathfrak{c}(\Omega,\mathcal{S},m)$ a bounded constant. Using the fact that $\vartheta<\mu\theta$, given in Lemma 2.2, we conclude from the above calculations that

[TABLE]

with $\widetilde{\mathfrak{c}}$ a bounded constant that depends on $\Omega$, $\phi$, $\mathcal{S}$, and $m$.

Since $j^{\prime}$ is taken arbitrarily, we conclude that

[TABLE]

By the symmetry in our calculations between the two groups of point sources, we see that we could replace the left hand side of the above inequality with the Hausdorff distance between the two groups of projected points $\{\mathcal{S}(\mathbf{x}_{1,j})\}_{j=1}^{m}$ and $\{\mathcal{S}(\mathbf{x}_{2,j})\}_{j=1}^{m}$:

[TABLE]

The stability result (11) then follows from this fact and the Hall theorem [Cameron-Book94, ElEl-CRASP12], which states that there exists a permutation $\pi$ acting on $\{1,2,\cdots,m\}$, that is a renumbering of the points, such that the Hausdorff distance can be realized by $\rho_{\mathbf{x}}$.

The next step is to establish the stability for the strengths of point sources. For an integer $1\leq j^{\prime}\leq m$, we introduce the function

[TABLE]

Then $\varphi_{j^{\prime}}$ solves the equation

[TABLE]

and satisfies

[TABLE]

Following the same procedure as before, we multiply the equation for $w$ by $\varphi_{j^{\prime}}$ and the equation for $\varphi_{j^{\prime}}$ by $w$, take the difference of the results, and apply Green’s identity to obtain,

[TABLE]

By the Cauchy-Schwartz inequality, we have

[TABLE]

where

[TABLE]

with $\vartheta$ and $\theta$ defined as before, and $\mathfrak{c}=\mathfrak{c}(\Omega,\mathcal{S},m)$ a bounded constant. We therefore have

[TABLE]

We now verify, using the assumption that $\sigma>\bar{\rho}_{\mathbf{x}}$, that,

[TABLE]

This allows us to conclude that, for some constants $\mathfrak{c}^{\prime}$ and $\mathfrak{c}^{\prime\prime}$ (for instance, one could take $\mathfrak{c}^{\prime}=\theta^{-1}\sigma^{1-m}(\sigma-\rho_{\mathbf{x}})^{1-m}$), we have

[TABLE]

The stability bound (12) then follows from  (18) and (19). ∎

The stability of reconstructing the locations of the point sources is Hölder type with exponent $\frac{1}{m}$, $m$ being the number of point sources included. The stability deteriorates quickly when $m$ increases. Therefore, we could only hope to reconstruct stably a very small number of point sources in practice.

The conditional stability of the reconstructing the strengths of the point sources contains two parts. The second part is from the Cauchy data and is Lipschitz type. The constant in front of it, however, depends on $m$. When ${\rm diam}_{\mathcal{S}}(\Omega)$ is large, this constant blows up quickly with $m$, another indication that one can not hope to stably reconstruct a large number of point sources. The first part comes from error in the determination of the locations of the point sources. If the locations are reconstructed perfectly, this term disappear. If, on the other hand, there is a large error in the reconstructing of the locations, the error in the reconstruction of the strengths is also large.

2.2 Stability with respect to media changes

Here we study the stability of the reconstruction of point sources with respect to smooth media changes. We assume that the measured data are collected with a medium $n_{1}$ that we do not know exactly. We then reconstruct the point sources pretending that the medium in which the data were collected is $n_{2}$. We show that the reconstructions in $n_{2}$ is not too different from the reconstructions in $n_{1}$ if $n_{2}$ is not too different from $n_{1}$, in appropriate sense.

Theorem 2.6**.**

Let $q_{1}$ and $q_{2}$ be two sources of the form (4) $(m:=m_{1}=m_{2})$, reconstructed for two media with refractive index $n_{1}$ and $n_{2}$ respectively, using Cauchy data $(f,g)$. Under the assumptions in (A) for $(q_{\ell},n_{\ell})$ ($\ell=1,2$), there exists a permutation $\pi$ acts on $\{1,2\dots,m\}$ such that

[TABLE]

$\sigma=\min_{i\neq j}|\mathcal{S}(\mathbf{x}_{2,i})-\mathcal{S}(\mathbf{x}_{2,j})|$, $\mu$ is from Lemma 2.2, and $\mathfrak{c}=\mathfrak{c}(\Omega,m,\overline{\lambda},g)$ is a bounded constant.

Proof.

Let $u_{\ell}$ $(\ell=1,2)$ be the solution to the Helmholtz equation (1) with the source and refractive index pair $(n_{\ell},q_{\ell})$. We define $w:=u_{1}-u_{2}$, $\delta n=n_{1}-n_{2}$, and $\delta q=q_{1}-q_{2}$. Then $w$ solves

[TABLE]

with boundary data $(w,\partial_{\nu}w):=(0,0)$.

Let $\phi$ be the function defined in Lemma 2.2 for the medium with refractive index $n_{1}$. Let $1\leq j^{\prime}\leq m$ be an integer. We use the function

[TABLE]

This function now solves the equation

[TABLE]

and satisfies

[TABLE]

Multiplying the equation for $w$ by $\psi_{j^{\prime}}$ and the equation for $\psi_{j^{\prime}}$ by $w$, taking the difference of the results, and applying Green’s identity, we have

[TABLE]

By the Cauchy-Schwartz inequality, we have

[TABLE]

where $K:=\operatorname{supp}(n_{2}-n_{1})\subset\subset\Omega$, $\|\psi_{j^{\prime}}\|_{L^{\infty}(\Omega)}\leq\|\phi\|_{L^{\infty}(\Omega)}({\rm diam}_{\mathcal{S}}(\Omega))^{2m-1}$.

Let $G$ be the fundamental solution for $\Delta+k^{2}(1+n_{2})$ with homogeneous Neumann boundary condition. We have the following representation for $u_{2}$:

[TABLE]

This allows us to conclude that $\|u_{2}\|_{L^{2}(K)}\leq\mathfrak{c}(1+m\overline{\lambda})$ for some constant $\mathfrak{c}=\mathfrak{c}(\Omega,K,g)$.

The rest of the proof is identical to that of Theorem 2.5. We first verify that

[TABLE]

where $\rho_{j^{\prime}}:=\displaystyle\min_{1\leq j\leq m}{\rm dist}_{\mathcal{S}}(\mathbf{x}_{2,j^{\prime}},\mathbf{x}_{1,j})$ is defined the same way as before.

Therefore, we have obtained

[TABLE]

The stability bounded (20) then follows from symmetry argument and the Hall theorem [Cameron-Book94, ElEl-CRASP12]. ∎

This result shows that the reconstruction of the locations of the point sources, up to a permutation, is relatively robust against changes in the underlying medium. However, the stability again deteriorates fast when the number of point sources increases.

Note that the above result is based on the assumption that we know exactly the number of point sources inside the medium. We do not have a general uniqueness result that allows us to determine the number of point sources from the measurement in this case. However, in some special cases, we can hope to reconstruct uniquely the point sources in a unknown medium that is has sufficiently simple structures, utilizing the fact that point sources are more singular compared to media variations.

Media with localized perturbations.

Let us consider the case where the medium $n_{2}$ is the homogeneous medium $n_{1}$ with a finite number of additional localized anomalies. That is, the refractive index $n_{2}$ is of the form:

[TABLE]

where $\Omega_{k}$ is the support of the $k$-th anomaly and $\tau_{k}$ is its strength. We assume that the point sources are away from the local anomalies of the medium, that is, ${\rm dist}(\mathbf{x}_{j},\Omega_{k})\geq\mathfrak{c}>0$ $\forall j,k$, for some $\mathfrak{c}$. Then we could follow the same proof in Lemma 2.1 to show that $w$ must vanish outside the support of $\delta n=n_{2}-n_{1}$. This allows us to uniquely (up to a permutation as before) determine the point sources that do NOT live in the support of the local anomalies of the medium. This provides a uniqueness argument for the numerical point source reconstructions in [BuLaTaTs-AFR09].

3 Numerical experiments

We now perform some numerical simulations in the context of the theoretical study in the previous section. We have two specific aims in mind: (i) when the underlying medium is known but heterogeneous, we want to see how well we can reconstruct point sources inside the medium; and (ii) when the underlying medium is not known, we want to see how the reconstruction of point sources are affected by the reconstruction of the medium.

We therefore intend to reconstruct the medium as well as the point sources inside the medium. We assume that we have access to multiple Cauchy data sets. We use a two-step reconstruction process. In the first step, we form differential data sets to eliminate the effect of the point sources and focus only on the medium. In the second step, we use a single data set to reconstruct the point sources, so that our simulations are in consistent with the theory developed in the previous section.

Media reconstruction.

Let $u$ be the solution to the Helmholtz equation (1), and $\widetilde{u}$ be the solution of the same equation but with boundary source $f(\mathbf{x})+h(\mathbf{x})$. Then we check that $v:=\widetilde{u}-u$ solves the Helmholtz equation

[TABLE]

By changing the probe source $h(\mathbf{x})$, we could obtain data determined by the Dirichlet-to-Neumann operator

[TABLE]

These data allow us to reconstruct the refractive index $n(\mathbf{x})$ since that is the only unknown quantity in (23). This inverse problem has been studied extensively; see for instance [Isakov-Book06, NaUhWa-IP13] and references therein.

We perform the reconstruction by reformulate the inverse problem as a minimization problem. Let us assume that we have data generated from $J$ different probe sources $\{h_{j}\}_{j=1}^{J}$. We reconstruct $n(\mathbf{x})$ by minimizing the following mismatch functional:

[TABLE]

where $g^{*}_{j}$ is the measured differential data corresponding to the probe source $h_{j}$ and the parameter $\beta$ is the strength of the regularization term.

We solve this minimization problem with a quasi-Newton method [NoWr-Book06, ReBaHi-SIAM06] where we use the adjoint state method to calculate the gradient of the objective functional with respect to the refractive index. Let $w_{j}$ $(1\leq j\leq J)$ be the solution to adjoint equation

[TABLE]

We can then show that the Fréchet derivative of $\Phi$ with respect to $n$ in direction $\delta n$ is given as

[TABLE]

In the minimization process, we solve the forward and adjoint Helmholtz problems (23) and (26) with a standard $P_{1}$ finite element solver.

Source reconstruction.

Once the refractive index $n(\mathbf{x})$ is reconstructed, we can reconstruct the unknown point sources, encoded in $q(\mathbf{x})$ in the Helmholtz equation (1), from observed boundary data. We do this again with a minimization strategy. More precisely, we minimize the functional

[TABLE]

over the locations and strengths of the point sources.

To avoid dealing with the singularity of the solution $u$ due to the point sources, we explicitly factorize out the singular part of $u$ as follows. Let $G(\mathbf{x};\mathbf{y})$ be the fundamental solution of the homogeneous Helmholtz operator in the whole space, that is,

[TABLE]

We represent the solution of (1) through the integral equation

[TABLE]

Let $\widehat{u}(\mathbf{x}):=u(\mathbf{x})+\sum_{j=1}^{m}\lambda_{j}G(\mathbf{x}_{j};\mathbf{x})$, then $\widehat{u}$ solves the integral equation

[TABLE]

where the source term

[TABLE]

To find the solution $u$, we solve for $\widehat{u}$ using (30) and then form $u=\widehat{u}-\sum_{j=1}^{m}\lambda_{j}G(\mathbf{x}_{j};\mathbf{x})$.

To evaluate the gradient of the objective function, we introduce the adjoint problem

[TABLE]

We can then show that the gradient of $\Psi$ with respect to a parameter the strength $\lambda_{k}$ and location $\mathbf{x}_{k}$ are given respectively as

[TABLE]

The numerical simulations we present below are all done in a two-dimensional domain for simplicity. The best way to make this consistent with the theory in the previous section, which are constructed in dimension three, is to view the the simulations as simplifications of three-dimensional ones for which the refractive index and the illumination sources are invariant in the $z$-direction. We set the domain $\Omega=[0,1]^{2}$ and the wave number $k=8$ in our experiments. We collect $J=6$ differential data sets generated from sources $f$ and $\{h_{j}\}_{j=1}^{J}$ to reconstruct the refractive index. To avoid the inverse crime, the synthetic measurements are generated on a fine grid while the inversion is fulfilled on another coarse grid. Moreover, we pollute our synthetic data with multiplicative random noise by perform the operation: $g^{*}_{j}\to g^{*}_{j}(1+\tau\mathcal{U}([-1,1]))$ with $\mathcal{U}[-1,1]$ the uniformly distributed random variable in $[-1,1]$ and $\tau$ the level of noise that we will specify later. The algorithms are implemented in the $\verb|MATLAB|$ software with the source codes deposited at $\verb|github|$ 111The github repository for our source codes is at https://github.com/lowrank/ips..

We performed simulations on several different media. Here we present results on two typical ones that have refractive indices respectively.

[TABLE]

where $\mathbf{y}_{1}=(0.25,0.25)$, $\mathbf{y}_{2}=(0.75,0.75)$, $R=0.25$ and $\chi_{D_{R}(\mathbf{y}_{k})}$ is the characteristic function of the disk of radius $R$ centered at $\mathbf{y}_{k}$, and

[TABLE]

where $\chi_{Rec}$ is the characteristic function of the rectangle $Rec=(0.5\ 0.75)\times(0.25\ 0.75)$; see Figure 3 and Figure LABEL:FIG:Unknown_Media_2 respectively for the plots of these refractive indices.

Experiment 1 [Recovery in Known Environments].

In the first set of numerical experiments, we perform reconstructions of point sources in heterogeneous media with known refractive indices. In Figure 1 and Figure 2, we show reconstructions of the locations of the point sources in the media (33) and (34) respectively. The true and reconstructed strengths are summarized in the first two rows of Table 3.