Inverse Problem of Diffraction by an Inhomogeneous Solid with a Piecewise Hoelder Refractive Index

TL;DR

This paper presents a two-step method for reconstructing an unknown piecewise-Hölder refractive index of an inhomogeneous solid from boundary measurements, proving uniqueness and verifying effectiveness through numerical tests.

Contribution

It introduces a novel two-step approach for inverse diffraction problems with piecewise-Hölder refractive indices, including proof of uniqueness and numerical validation.

Findings

Uniqueness of the current solution in piecewise-constant class.

Effective numerical reconstruction of the refractive index.

Validation of the method through test problems with accurate results.

Abstract

The problem of reconstruction of an unknown refractive index $k(x)$ of an inhomogeneous solid $P$ is considered. The refractive index is assumed to be a piecewise-H\"{o}lder function The original boundary value problem for the Helmholtz equation is reduced to the integral Lippman-Schwinger equation. The incident wave is defeined by a point source located outside $P.$ The solution of the inverse problem is obtained in two steps. First, the "current" $ J = (k^2 - k_0^2)u$ is determined in the inhomogeneity region. Second, the desired function $k (x) $ is expressed via the current $ J (x) $ and the incident wave $u_0.$ The uniqueness of the solution $J$ to the first-kind integral equation is proved in the class of piecewise-constant functions. The two-step method was verified by solving a test problem with a given refractive index. The comparison between the approximate solutions and the exact one approved the efficiency of the proposed method.

Equations

$$P=\{x=(x_1,x_2,x_3):a_1<x_1<b_1,a_2<x_2<b_2,a_3<x_3<b_3\}$$

$$x_{1,i_1}=a_1+\frac{b_1-a_1}{n}{i}_1,x_{2,i_2}=a_2+\frac{b_2-a_2}{n}{i}_2,x_{3,i_3}=a_3+\frac{b_3-a_3}{n}{i}_3,(0\le i_k\le n)$$

$${\Pi }_{i_1{i}_2{i}_3}=\{x:x_{k,i_k}<x_k<x_{k,k+1}\},0\le i_k\le n-1.$$

$${\chi }_{i_1{i}_2{i}_3}(x)=\{\begin{array}{ll}1,& x\in {\Pi }_{i_1{i}_2{i}_3},\\ 0,& x\notin {\Pi }_{i_1{i}_2{i}_3}.\end{array}$$

$$k(x)=\{k_{i_1{i}_2{i}_3}(x),x\in {\Pi }_{i_1{i}_2{i}_3},$$

$$k_{i_1{i}_2{i}_3}\in C^{0,\alpha }({\Pi }_{i_1{i}_2{i}_3}).$$

$$k(x)=\sum _I{k}_I(x){\chi }_I(x).$$

$${\Pi }_I^{\prime }={\Pi }_I\setminus E_P,P^{\prime }=P\setminus E_P.$$

$$U_0(x,t)=u_0{e}^{-i\omega t},U_s(x,t)=u_s{e}^{-i\omega t},U(x,t)=U_0(x,t)+U_s(x,t),$$

$$u_0(x)=\frac{e^{i{k}_0|x-x_0|}}{4\pi |x-x_0|},x_0\notin \overline{P}.$$

$$(△+k_0^2)u_0(x)=-\delta (x-x_0)$$

$$({\mathcal{P}}_1)\{\begin{array}{l}(△+k_I^2(x))u(x)=0,x\in {\Pi }_I;(△+k_0^2(x))u(x)=-\delta (x-x_0),x\in {\mathbb{R}}^3\setminus (\overline{P}\cup \{x_0\});\\ {[u]|}_{\partial {\Pi }_I}=0,{\left[\frac{\partial u}{\partial \mathbf{n}}\right]|}_{\partial {\Pi }_I^{\prime }}=0;\\ u\in H_{loc}^1({\mathbb{R}}^3\setminus \{x_0\});\\ \frac{\partial u_s}{\partial r}=i{k}_0{u}_s+o\left(\frac{1}{r}\right),(\mathrm{Im}{k}_0=0);u_s(r)=O\left(\frac{1}{r^2}\right),(\mathrm{Im}{k}_0>0).\end{array}$$

$$u\in C^1({\mathbb{R}}^3\setminus \{x_0\})\bigcap _I{C}^2({\Pi }_I)\bigcap C^2({\mathbb{R}}^3\setminus (\overline{P}\cup \{x_0\})),$$

$$\Delta u(x)+k_0^{2}u(x)=(k_0^2-k_I^2(x))u(x),x\in {\Pi }_I.$$

$$\Delta u(x)+k_0^{2}u(x)=-\delta (x-x_0),x\in {\Pi }_0,$$

$$\displaystyle\int\limits_{\partial\Pi_{I}}\left(\frac{\partial u(y)}{\partial\mathbf{n}}G(x,y)-\frac{\partial G(x,y)}{\partial\mathbf{n}}u(y)\right)ds_{y}=\int\limits_{\Pi_{I}}\Bigl{(}\Delta u(y)G(x,y)-\Delta G(x,y)u(y)\Bigr{)}dy=$$

$$\displaystyle=\int\limits_{\Pi_{I}}\Bigl{(}-k_{I}^{2}(y)u(y)G(x,y)+k_{0}^{2}G(x,y)u(y)+\delta(x-y)u(y)\Bigr{)}dy=$$

$$=u(x)-\underset{{\Pi }_I}{\int }(k_I^2(y)-k_0^2)G(x,y)u(y)dy,x\in {\Pi }_I;$$

$$\displaystyle\int\limits_{\partial\Pi_{J}}\left(\frac{\partial u(y)}{\partial\mathbf{n}}G(x,y)-\frac{\partial G(x,y)}{\partial\mathbf{n}}u(y)\right)ds_{y}=\int\limits_{\Pi_{J}}\Bigl{(}\Delta u(y)G(x,y)-\Delta G(x,y)u(y)\Bigr{)}dy=$$

$$\displaystyle=\int\limits_{\Pi_{J}}\Bigl{(}-k_{J}^{2}(y)u(y)G(x,y)+k_{0}^{2}G(x,y)u(y)\Bigr{)}dy=$$

$$=-\underset{{\Pi }_J}{\int }(k_J^2(y)-k_0^2)G(x,y)u(y)dy,x\in {\Pi }_I(J\ne I);$$

$$\displaystyle\int\limits_{\partial P\cup S}\left(\frac{\partial u(y)}{\partial\mathbf{n}}G(x,y)-\frac{\partial G(x,y)}{\partial\mathbf{n}}u(y)\right)ds=\int\limits_{B\setminus P}\Bigl{(}\Delta u(y)G(x,y)-\Delta G(x,y)u(y)\Bigr{)}dx=$$

$$\displaystyle=\int\limits_{B\setminus P}\Bigl{(}-k_{0}^{2}u(y)G(x,y)-\delta(y-x_{0})G(x,y)+k_{0}^{2}G(x,y)u(y)\Bigr{)}dy=-G(x,x_{0}).$$

$$\underset{S}{\int }\left(\frac{\partial u(y)}{\partial \mathbf{n}}G(x,y)-\frac{\partial G(x,y)}{\partial \mathbf{n}}u(y)\right)ds=u(x)$$

$$-G(x,x_0),x\in {\Pi }_I.$$

$$u(x)-\sum _J\underset{{\Pi }_J}{\int }(k_J^2(y)-k_0^2)G(x,y)u(y)dy=G(x,x_0),x\in {\Pi }_I.$$

$$u(x)-\underset{P}{\int }(k^2(y)-k_0^2)G(x,y)u(y)dy=u_0(x),x\in P$$

$$u(x)=u_0(x)+\underset{P}{\int }(k^2(y)-k_0^2)G(x,y)u(y)dy,x\in {\mathbb{R}}^3\setminus (P\cup \{x_0\}).$$

$$u(x)-\underset{{\Pi }_I}{\int }(k_I^2(y)-k_0^2)G(x,y)u(y)dy=\sum _{J\ne I}\underset{{\Pi }_J}{\int }(k_J^2(y)-k_0^2)G(x,y)u(y)dy+u_0(x),x\in {\Pi }_I.$$

$$u(x)=\underset{{\Pi }_I}{\int }(k_I^2(y)-k_0^2)G(x,y)u(y)dy+w(x)=v(x)+w(x),x\in {\Pi }_I,$$

$$c(y)=\{\begin{array}{ll}1,& y\in B=B_{d/4}(x_0),\\ 0,& y\in B={\Pi }_I\setminus B_{3d/4}(x_0).\end{array}$$

$$v(x)=v_1(x)+v_2(x)=\underset{{\Pi }_I}{\int }(k_I^2(y)-k_0^2)G(x,y)u(y)c(y)dy+\underset{{\Pi }_I\setminus B}{\int }(k_I^2(y)-k_0^2)G(x,y)u(y)(1-c(y))dy.$$

$$\begin{array}[]{ll}(\Delta+k_{e}^{2})v_{I}(x)=0,&x\notin P\;(I=-1,0),\\ (\Delta+k^{2}_{I}(x))v_{I}(x)=0,&x\in\Pi_{I},\end{array}$$

Full text

**Inverse Problem of Diffraction by an Inhomogeneous Solid with a Piecewise Hölder Refractive Index

**M.Yu. Medvedik, Yu.G. Smirnov and A.A. Tsupak

(Penza State University, Penza, Russia)

e-mail: smirnovyug@mail.ru; altsupak@yandex.ru

Introduction

Nowadays, when solving inverse problems of electrodynamics, one meets a problem of obtaining accurate solutions using a relatively small number of measurements. The present work is devoted to the solution of this very problem.

The urgency of the problem solved is due, on the one hand, to important applications, e.g., microwave tomography for early of breast cancer and reconstruction of the characteristics of samples of composite anisotropic materials.

On the other hand, methods and algorithms for solving inverse problems are important from the theoretical point of view in mathematical physics.

New numerical adaptive methods, parallel computational algorithms and supercomputer computations are among the most promising approaches to solving the inverse problems.

In particular, there are several approaches and many papers devoted to solving inverse problems of electrodynamics. Among theoretical works, in which the questions of the existence and uniqueness of solutions to the inverse problems are considered, we only mention the work [1], since it contains an extensive and up-to-date bibliography on the subject (63 references).

The most important methods for solving inverse problems are those that can be numerically implemented. These methods are mainly based on solving the hyperbolic systems of differential equations in the time domain using the finite differences or finite elements methods, with subsequent minimization of the corresponding functionals and Tikhonov regularization. Such methods and approaches are described in the monographs [2] - [8].

The volume singular integral equations method is an alternative approach to solving the inverse problems of electrodynamics. It is not widely used in electrodynamics (apparently, due to the greater complexity of numerical implementation in comparison with the finite element or finite difference methods. Nevertheless, it is successfully applied by both the authors of the paper and other researchers (e.g., [9]–[14]) to solving (direct) problems of diffraction of electromagnetic waves by dielectric solids.

The method of volume singular integral equations was applied, e.g., in [9] and [15], where the problems of objects reconstruction were solved using asymptotic data (far zone fields).

In the present paper, we use the method of integral equations for reconstruction of an unknown refractive index $k(x)$ (this function describes an inhomogeneous volume obstacle $P$ of a monochromatic wave). To solve the problem, we use a finite number of values of the given scalar field in the near zone, i.e., at several points in some region $D,$ lying outside the scatterer $P.$ Problems exactly like this one arise, e.g., in microwave tomography.

The article consists of two parts, the first of which is devoted to theoretical investigation of the problem and the description of the proposed method for its solution. The second part contains the description of computational experiments. Solving the inverse problem we pay considerable attention to the investigation of the direct diffraction problem. This is necessary for proving the equivalency between the differential and the integral formulations (the latter is used to solve the inverse problem) of the inverse problem of reconstructing the refractive index.

The first section of the article is devoted to the investigation of a direct problem diffraction of an external field by a volume body $P.$ The solid $P$ is characterized by a given refractive index.

First, we consider a boundary value problem for the Helmholtz equation in the classical formulation. Then, the original problem is reduced to the integral Lippman-Schwinger equation with respect to the unknown total field $u$ in the inhomogeneity region $P.$ In addition, we define $u$ outside $P$ using the integral representation. We assume that the incident wave is defeined by a point source located outside $P,$ whereas the refractive index is assumed to be a piecewise-Holder function.

The Lippmann-Schwinger equation is convenient to be studied in the $L_{2}(P)$ space for two reasons. First, in this case the operator of the equation is a Fredholm operator with index zero. Second, we are eager to use piecewise constant basis functions for numerical implementation of the method.

Investigating the integral equation in such a wide space, we show that, given a sufficiently smooth right-hand side, one obtains a smooth solution to the equation which, furthermore, represents a classical solution of the diffraction problem.

Thus, we show that the boundary value problem for the Helmholtz equation is equivalent to a system of integral equations, which is used to reconstruct the refractive index in the second section of the article.

We find the solution to the inverse problem using the two-step method (TSM). First, we determine the ”current” $J=(k^{2}-k_{0}^{2})u$ in the inhomogeneity region $P$ using the values of the total field given in the domain $D.$ For this end, we solve the integral equation of the first kind. It is shown in the paper that the solution to the integral equation is, in general, not unique. However, we prove that the solution is unique in the class of piecewise-constant functions. Such a choice of the solution space is seemingly valid for solving the inverse problem. Indeed, we are eager to use the collocation method with piecewise constant basis functions with compact support. In addition, the complete field $u$ corresponding to the piecewise-constant $J$ satisfies all the smoothness conditions that were formulated in the direct problem. That is due to the representation of the field $u$ via the volume potential with the density $J.$

In the second step of the proposed method, we express the desired function $k(x)$ via the current $J(x)$ and the incident wave $u_{0}$ using the Lippman-Schwinger equation in the domain $P.$

In the proposed method, we immediately define (and fixate) the computation grid. Therefore, classical theorems on the convergence of the method are not considered in the article. We believe that such an approach is valid in practice since it is usually clear what accuracy (and, consequently, what grid) is needed. E.g., it is sufficient to consider grids with 5mm step in the method of microwave tomography for the breast cancer early diagnosis (thus, any grid refinement is no use).

The main advantage of the proposed method is the reduction of the original boundary value problem to the solution of the  emph linear integral equation with the subsequent explicit calculation of the desired function. Thus, we are spared, first, the need to solve any nonlinear equations. Secondly, we don t apply iterative methods that require good choice of an initial approximation.

One of the main results of the present paper is the proof of the uniqueness theorem for the solution to the integral equation of the first kind in the class of piecewise-constant functions. In addition, we show that even in the class of analytic functions the solution is not unique.

The main difficulty of the numerical method implementation is in solving the integral equation of the first kind with a smooth kernel. Note, that we apply collocation method. The corresponding matrices are ill-conditioned (albeit nondegenerate), which implies a certain instability with respect to the right-hand side of the equation. As a result, given distorted original data (i.e., the given values of the field in the domain $D$), one obtains ”false” inhomogeneities in the region $P.$ In this paper, such inhomogeneities are called ”artifacts” of the reconstructed function $k(x).$ The computational experiments showed that the ”artifacts” emerged when the domain $D$ was moved away from the body $P.$

The authors propose two ways to eliminate the ”artifacts.”

The first method consists in the screening out ”extraneous” noise, i.e. in filtering the input data of the integral equation. When implementing this procedure, we assume that the ”true” noise should be a small (but sufficiently smooth) disturbance of the field. Data that does not satisfy this condition is discarded.

The second method consists in changing the location of the field source (another choice of the point $x_{0}$) and the receivers (these are the collocation points in the domain $D$). As a result, we determine the true inhomogeneities in the region $P$ and discard the false ones (the ”artifacts”). Note that the artifacts may change their location. In addition, the ”artifacts” are located in the $P$ region quite symmetrically which indirectly confirms their ”artificial” nature. Application of the rotation method allows to restore the inhomogeneity of the body with sufficient accuracy.

To refine the solution, we apply the following adaptive method: in the region of the inhomogeneity found, we define new (and more dense) grid, then carry out additional measurements. As a result, a more accurate value of the refractive index is obtained.

The authors solved a series of problems, analyzed the solutions obtained and compared them with the exact solutions. A description of one of the numerical experiments is presented in this article.

1 The Direct Scattering Problem

Prior to the describing the statement of the inverse problem as well as the method for its solving, we shall investigate the direct diffraction problem.

1.1 The Boundary Value Problem of Diffraction by a Solid with a Piecewise Hölder Refractive Index

Consider an isotropic inhomogeneous rectangle parallelepiped

[TABLE]

located in the homogeneous space $\mathbb{R}^{3}.$

Define the uniform mesh in the domain $\overline{P}$ as follows

[TABLE]

and introduce sub-domains $\Pi_{i_{1}i_{2}i_{3}}:$

[TABLE]

We also define a set of piecewise constant functions $\chi_{i_{1}i_{2}i_{3}}$ (indicator functions):

[TABLE]

Henceforward we assume that inhomogeneity of the domain $P$ is described by a piecewise continuous function $k(x)=n(x)k_{0}$ such that

[TABLE]

where all functions $k_{i_{1}i_{2}i_{3}}(x)$ are Hölder continuous

[TABLE]

Note that at the points of the parallelepipedal faces $\partial\Pi_{i_{1}i_{2}i_{3}}$ the function $k(x)$ can be defined via one-sided limits from either side of the face.

By introducing the multi-indices $I=(i_{1}i_{2}i_{3})$ one can now define the function $k(x)$ for any point $x\in P$ by the following equality:

[TABLE]

The lossless medium outside the solid is characterized by a given positive wavenumber $k_{0}>0.$

Define $E_{P}$ as the union of all edges of the parallelepipeds $\Pi_{I}$ and give the following notation:

[TABLE]

The incident wave (the source field) as well as the scattered and the total fields are considered to depend on time harmonically:

[TABLE]

Thus it is sufficient to formulate the scattering problem for the scalar complex amplitude $u(x)$ of the total field.

We consider the incident field of a point source setting

[TABLE]

The field represents a solution to the Helmholtz equation

[TABLE]

that satisfies Sommerfeld radiation conditions.

The direct scattering problem in the rigorous mathematical statement is to find a solution $u(x)$ to the following boundary value problem:

[TABLE]

Definition 1**.**

Any solution to the problem $({\cal P}_{1})$ that satisfies the conditions

[TABLE]

of continuity is a quasiclassical solution to the direct scattering problem.

1.2 Lippman-Schwinger Integral Equation. Smoothness of Solutions to the Integral Equation.

Now we are eager to reduce the problem $({\cal P}_{1})$ to the Lippman-Schwinger integral equation.

Rewrite the Helmholtz equation in sub-domains $\Pi_{I}$ as follows

[TABLE]

In the bounded region $\Pi_{0}=B\setminus\overline{P}$ one obtains

[TABLE]

where $B\supset P$ is a sufficiently large ball centered at zero (denote also $S=\partial B$) and $G(x,y)=\frac{\exp(ik_{0}|x-y|)}{4\pi|x-y|}$ is the Green function of the Helmholtz equation. Applying the second Green formula one derives

[TABLE]

[TABLE]

[TABLE]

Add equalities (8)-(10) and take into account the transmission conditions:

[TABLE]

Passing in (11) to the limit as the radius of the ball $B$ tends to infinity one gets

[TABLE]

The last equation can be rewritten as follows

[TABLE]

since the definition of the function $k(x).$ Consider also the integral representation of the total field in the outside of the solid $P:$

[TABLE]

Definition 2**.**

The integral statement of the direct diffraction problem is understood as the system $({\cal P}_{2})$ consisting of equation (13) in the domain $P$ and representation (14) outside it.

The operator in equation (13) is denoted by ${\cal I}-{\cal A}$ and is treated as a mapping in the $L_{2}(P)$ space.

First, let us show that any solution $u(x)$ of the problem $({\cal P}_{2})$ satisfies the smoothness conditions, formulated in the quasiclassical statement of the problem.

Theorem 1**.**

Let equation (13) have a solution $u\in L_{2}(P).$ Then, the smoothness conditions (5) are satisfied by the total field $u(x),$ extended outside $P$ according to (14).

Proof.

From the definition of the incident wave in the considered statement of the problem it follows that $u_{0}\in C^{\infty}(\mathbb{R}^{3}\setminus\{x_{0}\}).$

Any solution is infinitely differentiable outside the solid since the smoothness of the integral operator kernel at each $x\notin\overline{P}.$

Now consider equation (13). For each multi-index $I$ one has

[TABLE]

The righthand side of the equation is infinitely differentiable in the open domain $\Pi_{I}$ since, for any $J\neq I,$ one has $G(x,y)\in C^{\infty}(\Pi_{I}\times\Pi_{J}).$

The inclusion $\in L_{2}(P)$ implies $u\in H^{2}(P)$ and, consequently, $u\in C^{\alpha}(\overline{P})$ for all $0<\alpha<1/2.$ Then ${\cal A}u\in C^{1}(\mathbb{R}^{3}\setminus\{x_{0}\})$ (see. [16]). From the latter follows the inclusion $u\in C^{1}(\mathbb{R}^{3}\setminus\{x_{0}\})$ which results also in the energy finiteness condition $u\in H^{1}_{\mathrm{loc}}(\mathbb{R}^{3}\setminus\{x_{0}\}).$

It suffices to prove that $u\in C^{2}(\Pi_{I})$ for any $I.$ Write the following equality for $u$ according to equation (15):

[TABLE]

where $w\in C^{\infty}(\Pi_{I}).$

Let $x_{0}\in\Pi_{I}$ be an arbitrary inner point of the $I$-th sub-domain such that $d=dist(x_{0},\partial\Pi_{I})>0.$ Introduce the cut-off function $c\in C^{\infty}(\Pi_{I}):$

[TABLE]

Represent $v$ in the form given below:

[TABLE]

Since the smoothness of the kernel in the second term, one obtains $v_{2}\in C^{\infty}(B).$

Note that $u\in C^{0,\alpha}(\Pi_{I})$ and $c\in C_{0}^{\infty}(\Pi_{I}),$ $(k^{2}_{I}(y)-k_{0}^{2})G(x,y)u(y)c(y)\in C_{0}^{0,\alpha}(\mathbb{R}^{3}).$ From the latter, using the properties of the volume potential (for details, see [9] on page 212), follows the inclusion $v_{1}\in C^{2,\alpha}(\mathbb{R}^{3}).$

Thus, a solution $u(x)$ is twice differentiable in a vicinity of each point $x_{0}\in\Pi_{I},$ i.e. $u\in C^{2}(\Pi_{I}).$

∎

1.3 The Equivalency Theorem. Uniqueness of a Solution to the Diffraction Problem

Let us formulate and prove two important results of investigation of the direct diffraction problem, which are the theorem on equivalency between the differential in the integral formulations of the problem, and the theorem on uniqueness of its quasiclassical solution.

Theorem 2**.**

The problems $({\cal P}_{1})$ and $({\cal P}_{2})$ are equivalent. More precisely, if $u(x)$ is a quasiclassical solution to the problem $({\cal P}_{1})$ then $u$ satisfies equation (13) and representation equality (14). Vise versa, for any solution $u\in L_{2}(P)$ to the integral equation (13), the total field $u(x),$ extended to $\mathbb{R}^{3}\setminus\{x_{0}\}$ by formula (14), is a quasiclassical solution to the problem $({\cal P}_{1})$.

Proof.

The former part of the theorem follows from the derivation of the integral equation.

Let $u$ be a solution to equation (13) with $u_{0}\in C^{\infty}(\mathbb{R}^{3}\setminus\{x_{0}\})$.

The definition of the term $u_{1}$ via a volume potential together with smoothness of the term $u_{0}$ in $\mathbb{R}^{3}\setminus\{x_{0}\}$ imply that $u$ is a solution to the Helmholtz equation in the domains $\Pi_{I}$ and $\mathbb{R}^{3}\setminus(\{x_{0}\}\cup\overline{P}).$

The scattered field $u_{s}(x)=\int\limits_{P}(k^{2}(y)-k_{0}^{2})G(x,y)u(y)dy$ satisfies the radiation condition, whereas the transmission conditions are fulfilled since the inclusion $u\in C^{1}(\mathbb{R}^{3}\setminus\{x_{0}\})$ shown in Theorem 1. Note that the equality $\frac{\partial u}{\partial\mathbf{n}}=0$ can not be written on the edges of the sub-domains $\Pi_{I}.$ ∎

Theorem 3**.**

For any $\mathop{\mathrm{Im}}\nolimits k(x)\geq 0$ the problem $({\cal P}_{1})$ has at most one quasiclassical solution.

Proof.

Let us show that the corresponding homogeneous boundary value problem (with $u_{0}\equiv 0$ in $\mathbb{R}^{3}$) formulated for the scattered field $u_{s}\equiv u$ has only the trivial solution $u_{s}\equiv 0.$

1. Consider a sufficiently large ball $B\supset\overline{P}$ of a radius $R.$ Introduce regions $\Pi_{0}:=B\setminus(\overline{P})$ with the boundary $\partial\Pi_{0}=\partial B\cup\partial P,$ and $\Pi_{-1}:={\overline{B}}^{c}.$

We shall reduce the original problem for the scattered field $u_{s}$ to a transmission problem in the domains $\Pi_{i}.$ For this purpose we denote the restriction of $u_{s}(x)$ to the closed subsets $\overline{\Pi}_{I}$ by $v_{I}(x).$ The functions $v_{I}$ satisfy the Helmholtz equation in the corresponding domains:

[TABLE]

as well as the transmission conditions at the boundaries of the adjacent domains $\Pi_{I},\Pi_{J}$:

[TABLE]

The radiation condition is now formulated for the function $v_{-1}:$

[TABLE]

The first Green formula, applied to the functions $\overline{v}_{I},v_{I}$ in the bounded domains $\Pi_{0}$ and $\Pi_{I}\subset P,$ as well as the Helmholtz equation yield in the following relation:

[TABLE]

Add equalities (20), using the transmission conditions (18):

[TABLE]

Consider the imaginary part of the latter relation and take into account the radiation conditions:

[TABLE]

Application of the Rellich lemma (see [17] on p.88) results in equality $v_{-1}(x)\equiv 0$ for each point $x\in\Pi_{-1}.$

3. Show now that the relation $v_{-1}(x)\equiv 0$ holds at the points of the inhomogeneity domain.

Consider an arbitrary ¡¡external¿¿ parallelepiped $\Pi_{I}$ such that $\partial\Pi_{I}\cap\partial P=S\neq\varnothing.$ Represent the solution $u$ as follows (see also the previous subsection):

[TABLE]

It can be shown that the function $u(x)$ is infinitely differentiable in a sufficiently small vicinity $U$ of an arbitrary point $x_{0}\in S.$ Consider the following representation of $w(x):$

[TABLE]

where $c(y)\in C^{0,\alpha}_{0}(U)$ is a compactly supported cut-off function such that $c(y)\equiv 1$ in $B_{r}(x_{0})=U^{\prime}\subset U.$ Then the smoothness of the integral operator’s kernel implies $w_{2}\in C^{\infty}(U^{\prime}).$ Further, the inclusion $w_{1}\in C^{2}(\mathbb{R}^{3})$ is valid since the term $w_{1}$ is a Newtonian potential with the compactly supported smooth density $(k^{2}(y)-k_{0})u(y)c(y)\in C_{0}^{0,\alpha}(\mathbb{R}^{3})$ (see [9] on p.207).

Thus, the function $u\in C^{2}(U^{\prime})$ is a solution to the Helmholtz equation such that $u\equiv 0$ in the sub-domain $U^{\prime}\setminus P.$ From the unique continuity principle ([9], p.212) it follows now that $u\equiv 0$ in the domain $U^{\prime}$ and, consequently in the entire parallelepiped $\Pi_{I}.$

Similarly, one can consider all sub-domains $\Pi_{I},$ repeat the above arguments, and conclude that $u\equiv 0$ in $P.$

If the condition $\mathop{\mathrm{Im}}\nolimits k>0$ holds in the entire space, then, from the second relation in (19) it follows that $u(x)=O(R^{-2})$ an the sphere $\partial B.$ Consequently, the left-hand side of equality (21) vanishes as $R\to+\infty.$ As a results, one obtains

[TABLE]

For the imaginary part of (24) holds the following relation:

[TABLE]

as $R\to+\infty.$ Both terms in the latter expression are of the same sign due to the properties of the medium. From that we conclude that $v_{I}(x)\equiv 0,\;x\in\Pi_{I}.$

If $k_{0}>0$ outside the inhomogeneity domain and $\;\mathop{\mathrm{Re}}\nolimits k\!\cdot\!\mathop{\mathrm{Im}}\nolimits k(x)>0$ inside it, then we similarly deduce that $u\equiv 0$ in $P:$ using the Rellich lemma we first get $v_{-1}\equiv 0,$ and then, as in the item 3. of the current proof, we obtain that the solution is trivial inside the solid $P.$ ∎

The next statement on results from Theorems 2 and 3:

Theorem 4**.**

The operator

[TABLE]

is continuously invertible

Proof.

For any $u\in L_{2}(P)$ one has ${\cal A}u\in H^{2}(P).$ From the latter it follows that ${\cal A}:L_{2}(P)\to L_{2}(P)$ is a compact operator.

Let $u_{0}\equiv 0$ in $\mathbb{R}^{3}.$ Then the boundary value problem $({\cal P}_{1})$ has only the trivial solution (see Theorem 3). Thus, due to the equivalency between $({\cal P}_{1})$ and $({\cal P}_{2}),$ $u=0$ is the only solution to the integral equation $({\cal I}-{\cal A})u=0.$

Thus, $({\cal I}-{\cal A}):L_{2}(P)\to L_{2}(P)$ is an injective Fredholm operator with index zero. ∎

2 The Inverse Problem of Reconstructing the Refractive Index

2.1 Statement of the Inverse Problem

Consider in $\mathbb{R}^{3}$ an inhomogeneous parallelepiped $P$ characterized by an unknown refractive index $n(x).$ Assume, as in the statement of the direct diffraction problem, that $n(x)$ and $k(x)=n(x)k_{0}$ are piecewise-Hölder functions in the domain $P$ with a given mesh and a set of sub-domains $\Pi_{I}.$

Introduce a bounded domain $D$ such that $\overline{D}\cap\overline{P}=\varnothing,$ and assume that in the points $x\in D$ we are given the known values of the total field at a fixed frequency $\omega:$

[TABLE]

The monochromatic incident wave $U_{0}(x,t)=u_{0}(x)e^{-i\omega t}$ is defined according to (4), and the source of the field is located in an arbitrary point $x_{0}\notin\overline{P}\cup\overline{D}.$

In the proposed statement of the inverse diffraction problem we use the system $({\cal P}_{2})$ of integral equalities which represent the relation between the total field $u(x)$ and the function $k(x)$ (it is shown above that formulations $({\cal P}_{1})$ and $({\cal P}_{2})$ are equivalent).

We are eager to reconstruct the function $k(x)$ in the parallelepiped $P$ using measurements of the total field $u(x)$ at points of the bounded domain $D:$

[TABLE]

and taking into account the equation

[TABLE]

in the inhomogeneity domain $P.$

2.2 The Two-step Method for Solving the Inverse Diffraction Problem

In the domain $P,$ we introduce the function

[TABLE]

assuming that the condition $|k(x)|\geq\tilde{k}>k_{0}$ holds everywhere in $P.$ From the representation of the total field in the outside of the solid $P$ follows the equation

[TABLE]

for determination of $J(x),$ where as equation (13) can be rewritten as below:

[TABLE]

The idea of the proposed two-step method for reconstruction of the unknown coefficient $n(x)$ is as follows:

• •

Given the known values of the incident wave $u_{0}(x)$ and the total field $u(u)$ in the domain $D,$ we find, in the domain $P,$ the solution $J$ to equation (28).

• •

We reconstruct the function $k(x)$ inside $P$ using relation (29).

2.3 On Non-uniqueness of a Solution to the Integral Equation

Let us show that the homogeneous integral equation (28) has non-trivial solutions for any $k_{0}.$ For example we give below an argumentation for the case of a cubic solid $P=[-1;1]^{3}\subset R^{3}.$

Consider the function $\psi(x)=(1-x_{1}^{2})^{2}(1-x_{2}^{2})^{2}(1-x_{2}^{2})^{2}$ and introduce $J(x)=-(\triangle+k_{0}^{2})\psi.$ As $\psi$ satisfies the homogeneous boundary conditions \psi\bigr{|}_{\partial P}=\frac{\partial\psi}{\partial\mathbf{n}}\bigr{|}_{\partial P^{\prime}}=0, then the representation

[TABLE]

can be given. Introduce the potential $v(x)=\int\limits_{P}G(x,y)J(y)dy,\quad x\in\mathbb{R}^{3}.$ Then at any point $x\in P$ one obtains $v(x)=\psi(x).$ However, the relation $v\equiv 0$ holds outside the closed cube $\overline{P}:$

[TABLE]

The similar result can be obtained in the case of a domain $P$ of an arbitrary shape. For that purpose, define $\psi$ as an arbitrary smooth function with a compact support in $P$ so as to satisfy the conditions \psi\bigr{|}_{\partial P}=\frac{\partial\psi}{\partial\mathbf{n}}\bigr{|}_{\partial P^{\prime}}=0 and repeat the above analysis.

2.4 On Uniqueness of a Piecewise-constant Solution $J(x)$

It is shown above that the integral equation of the first kind has an infinite set of smooth solution.

We prove below that a unique solution $J$ can be obtained in the class of piecewise constant functions. Note that for approximate solving equation (28) we are going to apply the collocation method. That is why the choice of piecewise constant functions $J$ is a reasonable one. Thus, we find $J$ in the following way:

[TABLE]

where $J_{I}\in\mathbb{C}$ are the unknown coefficients, and $\chi_{I}(x)$ are the indicator functions of the sub-domains $\Pi_{I}$ (at points $x\in\Pi_{I}$ the function $J(x)$ can be defined by any constant value). Note that the class of piecewise constant solutions is sufficient for solving applied problems of physics, medical tomography etc. (see the introduction of the paper).

Below, we formulate and prove the theorem on uniqueness of a piecewise constant solution $J(x)$ to equation (29).

Theorem 5**.**

Consider a fixed set of $n^{3}$ rectangular sub-domains $\Pi_{I}$ in the inhomogeneity domain $P.$ Let

[TABLE]

If equation

[TABLE]

has a piecewise constant solution $J(x),$ then this solution is the unique one.

Proof.

1. Consider the homogeneous equation

[TABLE]

Introduce the volume potential

[TABLE]

Since $J(x)$ is a piecewise constant function then $v(x)\in C^{1}(\mathbb{R}^{3}).$ As a result, the transmission conditions

[TABLE]

hold on the boundaries of the parallelepipeds $\Pi_{I}.$

Moreover, the inclusions $v\in C^{2}(\Pi_{I})$ are valid. Consequently, the Helmholtz equation

[TABLE]

holds in the inner points $x\in\Pi_{I}$ (in the classical sense).

Outside $\overline{P},$ we obtain $(\triangle+k_{0}^{2})v(x)=0$ and $v\in C^{\infty}(\mathbb{R}^{3}\setminus\overline{P}).$

By the assumption of the theorem, the function $v$ is equal to zero in the domain $D\subset\mathbb{R}^{3}\setminus\overline{P}.$ Then, applying the unique continuation principle ([9], p.212), we derive that $v\equiv 0$ everywhere in $\mathbb{R}^{3}\setminus\overline{P}.$

From the inclusion $v(x)\in C^{1}(\mathbb{R}^{3})$ follows the relation

[TABLE]

Introduce the fundamental solution $\bar{G}(x,y)=\frac{e^{-ik_{0}|x-y|}}{4\pi|x-y|}$ to the Helmholtz equation. Apply the second Green formula to the functions $v,\bar{G}$ in the domains $\Pi_{I},$ take into account the homogeneous boundary conditions (37) and the transmission conditions on $\partial\Pi_{I}$ for the function $v:$

[TABLE]

From the latter follows the equality

[TABLE]

Subtracting (39) from equation (34) one deduces

[TABLE]

As the potential $v$ so the function $w$ satisfies the transmission conditions on $\partial P$ which implies

[TABLE]

In addition, the function $w\in C^{\infty}(\mathbb{R}^{3}\setminus\overline{P})$ satisfies outside $\overline{P}$ the equation $(\triangle+k_{0}^{2})w=0$ and the Sommerfeld radiation conditions. Consequently, $w\equiv 0$ in $\mathbb{R}^{3}\setminus\overline{P}.$

3. Thus, $w\equiv 0$ in $\mathbb{R}^{3}.$ Then for the Fourier transform ${\cal F}w(\xi)$ holds the similar relation ${\cal F}w\equiv 0$ everywhere in $\mathbb{R}^{3}.$

Introduce the grid parameters

[TABLE]

and the parallelepiped

[TABLE]

Then all other finite elements can be defined via shifts of the sub-domain $\Pi_{0}$ by appropriate vectors:

[TABLE]

The function $w(x)$ can now be represented as follows

[TABLE]

Evaluate the Fourier transform of $w$ taking in to account the latter relation:

[TABLE]

The relation ${\cal F}w\equiv 0$ is reduced to the equality

[TABLE]

on the centered sphere $S_{k_{0}}$ of the radius $k_{0}.$

By $S^{n-1}$ we denote a unit sphere in $\mathbb{R}^{n}$ ($n\geq 2$).

4. Let us shown that the functions $e^{ir_{I}\cdot\xi}$ are linearly independent on the sphere $S_{k_{0}}.$ To this end, we shall prove that the corresponding Gram matrix $\Gamma$ is nonsingular.

For an arbitrary matrix element $\Gamma_{II^{\prime}}$ one deduces

[TABLE]

In the above evaluation that the integrals $\int\limits_{S^{n-1}}f(\omega\cdot\xi)ds_{\xi}$ over the unit centered sphere $S^{n-1}$ do not depend on the variable $\omega\in S^{n-1}$ and, consequently, can be presented [18] as below:

[TABLE]

From (45) it follows that

[TABLE]

5. Represent $\Gamma$ via sum

[TABLE]

where $\tilde{I}$ is the unit matrix, and obtain the estimate

[TABLE]

The latter implies that the determinant of the Gram matrix is nonzero, since its diagonal element dominate at sufficiently large $n.$

Fix the row index $I^{\prime}=(0,0,0)$ and define $h=\min\{h_{1},h_{2},h_{3}\}$:

[TABLE]

The proof is complete. ∎

Since equation (29) represents the relation between $J,$ $k$ and $u,$ it now follows from theorem 5 that the solution $k(x)$ to the inverse diffraction problem, corresponding to the piecewise constant $J,$ is also unique.

Remark 1* (on existence of a solution).*

Let us assume that the righthand side of equation (32) is an element of a linear span of functions $\int_{\Pi_{I}}G(x,y)dy$ (for a given mesh on $p$). Then the operator of the lefthand side of equation (32) can be treated as mapping in finite-dimensional spaces. Such mapping are continuously invertible which results from theorem 5. In addition, the solutions to the inverse problem depend continuously on the given data.

2.5 Numerical Solution of the Inverse Problem

To solve equation (28) numerically, we apply the collocation method.

The current $J$ is sought in the form of a linear combination $\sum\limits_{j=1}^{N}c_{j}v_{j}(x)$ of piecewise constant basis functions. The collocation points are defined as follows: $r_{i}\in D$ ($i=1,\ldots,N.$)

Below we describe the conditions and results of computational experiments.

We consider the cubic inhomogeneity region $P$ with the edge length of $0.15m.$ The area of such a size is quite consistent with the objects studied, e.g., in breast cancer diagnosis.

There was carried out a series of computational experiments under various conditions. First, the array of receivers (points for measuring the total field) was defined on planes parallel to $xy-$ or $xz-$ planes. Second, the distance from the nearest receiver to the scatterer was varied (this parameter is denoted by $d_{r}$). Third, the calculations were made with or without artificially introduced errors (noise). The distance $d_{s}$ from the source of the incident wave $u_{0}$ to the body $P$ is fixed: $d_{s}=0.003m.$

All the figures below graphically represent the values of the real (under the letter (a)) and imaginary (under the letter (b)) parts of the function $k(x).$

We consider the following sample problem (SP): the inhomogeneity domain $P$ is characterized by a given complex-valued function $k(x).$ The graphical representation of the real and imaginary parts of the function $k(x)$ is in the figures 1(a) and 1(b), respectively.

Given the function $k(x),$ we first solve the direct diffraction problem. As a result, the field $u(x)$ and the ”current” $J(x)$ are determined in the region $P.$ These functions are then used for modeling the total field of the inverse problem.

Finally, the solution to the inverse is found using the two-step approximate method. The functions $\mathop{\mathrm{Re}}\nolimits k(x)$ and $\mathop{\mathrm{Im}}\nolimits k(x),$ shown in figures 2(a) and 2(b), correspond to the exact input data (no noise is added). This data is the values of the scattered field $u_{s}=u-u_{0}$ at the receivers points. Here $u_{0}$ is the given incident wave and $u$ is the solution to the direct problem of diffraction.

In the next experiment (see Figures 3(a) and 3(b)), the distance from the body to the receivers is increased up to $0.05m.$ In this case, the approximate solution significantly differs from the exact one. One can see numerous ”artifacts” in the figures.

In the third experiment (see Figures 4(a) and 4(b)), the distance $d_{r}$ is reduced to $0.005m,$ whereas the data is noisy. Here again, the solution is quite inaccurate and shows the ”artifacts”. In addition, there is also a loss of ”true” heterogeneities.

In the fourth experiment (see Figures 5(a) and 5(b)), the distance $d_{r}$ equals $0.005,m$ the input data is noisy, and the location of the receivers is changed.

In the final test (see Figures 6(a) and 6(b)), we changed the location of both the field source and the receivers. In this case, a more accurate solution was obtained. Despite some ”artifacts”, the values of the refractive index just slightly differ from the true ones.

Conclusion

The authors developed and theoretically justified a new method for solving the problem of reconstructing the refractive index in an inhomogeneous body using near-field data. The method involves solving a linear integral equation in the region of inhomogeneity.

The computational experiments confirm the efficiency of the developed algorithm. It is supposed to carry out experimental studies using the proposed method.