Numerical method for inverse scattering in two-layered background in near-field optics

Abstract

The inverse scattering problem in two-layered background that arises in near-field optics is discussed. The reconstruction of the scatterer from inhomogeneous medium deposited on a homogeneous substrate is presented, where the measured data only lies on a limited aperture. An error bound of the Born approximation is given. A moment method with truncated SVD is developed to handle the exponentially ill-posed problem with noisy measured data. Numerical experiments are presented to illustrate the resolution of the method, which are satisfying even under the perturbed data.

Keywords:

• inverse scattering
• two-layered background
• moment method
• truncated SVD
• near-field optics

1 Introduction

We can obtain images with subwavelength resolution via evanescent wave field in near-field optics.[1] There are three important modalities in near-field optics : near-field scanning optical microscopy (NSOM),[2] total internal reflection microscopy (TIRM)[3, 4] and photon scanning tunnelling microscopy (PSTM).[5, 6] The two principal genres of NSOM are illumination mode and collection mode. In illumination mode NSOM, a probe is scanned over the scatterer in the near zone and the scattering field is measured in the far zone of scatterer. In collection mode NSOM, the scatterer is illuminated by a source located in the far zone of the scatterer and the total field is collected in the near zone of the scatterer through the small aperture in the probe as the probe is scanned over the scatterer. In TIRM, the scatterer is illuminated by an evanescent wave that is generated by total internal reflection, and the scattering field is measured in the far zone of the scatterer. In the PSTM, the scatterer is illuminated by an evanescent wave and the scattering field is measured in the near zone.

Figure 1. The geometry of the model problem. The scatterer $q(\mathit{x})$ is in a compact support in $D$, the measurements are taken on $\mathrm{\Gamma }$ and $\mathrm{\Omega }$ is a compact set containing $D$ and $\mathrm{\Gamma }$.

In all of the three modalities, reconstruction of the scatterer from the measured data is desirable. In [1], within the framework of weak scattering, the authors give some numerical experiments to recover the scatterer. Bao and Li developed a recursive linearization method from an initial guess based on weak scattering. In [7], the model problem employs the first-order boundary condition, and requires the boundary measurements of scattering waves, but it only needs single-frequency incident wave. In [8], the radiation condition is employed, and it requires multi-frequency incident waves. The model in [8] uses a limited aperture data.

In this paper, we adopt the model with radiation condition, single-frequency incident wave and a limited aperture data. We will give an error estimate of the inverse problem. In [9], the authors deal with an optical diffusion problem, and give an inversion of the linear problem, but the computation is limited by the numerical precision, for its ill-posedness. Due to the ill-posedness of the inverse problem, we develop an effective numerical method: moment method with truncated SVD (singular value decomposition).

This paper is organized as follows. In Section 2, we give a mathematical model of the scattering problem. Next, based on the work of Schotland [10, 11], etc., the inverse problem is discussed. In Section 3, we discuss the regularization of the linearized problem. We present the moment method with truncated SVD, in order to handle the severely ill-posed problem with noisy measured data. Numerical examples are presented in Section 5 and some remarks are concluded in the final section.

2 Scattering and inverse scattering problem

2.1 Scattering problem

We describe the scattering model mathematically as follows. Assume that the materials are nonmagnetic in the problem considered in this paper. Throughout, we only consider the problem which can be reduced to a two-dimensional model in the case of transverse electric polarization. The index of refraction in the lower half-space has a constant value $n_0$. The index of refraction in the upper half-space varies within the domain of the sample but otherwise has a value of unity.[1, 8] That is, the background media consists of two layers with the index of refraction$$n(\mathit{x})=\left\{\begin{array}{cc}1,\hfill & \text{for}{x}_2>0,\hfill \\ n_0,\hfill & \text{for}{x}_2<0,\hfill \end{array}\right.$$where $\mathit{x}=(x_1,x_2)$ and $n_0>1$.

We denote the scatterer by $q(\mathit{x})$ within a compact support in $D\subset {\mathbb{R}}_{+}^2:=\{\mathit{x}=(x_1,x_2):x_2>0\}$, with $q(\mathit{x})/4\mathit{\pi }$ the susceptibility of the scatterer.[1] The incident wave comes from the lower half-space, which is a time-harmonic plane wave$$\begin{array}{c}\hfill u^i=e^{\mathrm{i}\mathit{\alpha }{x}_1+\mathrm{i}\mathit{\beta }{x}_2},\end{array}$$where$$\mathit{\alpha }=n_0{k}_{0}sin\mathit{\theta },\mathit{\beta }=n_0{k}_{0}cos\mathit{\theta },$$and $\mathit{\theta }\in (-\mathit{\pi }/2,\mathit{\pi }/2),$$k_0$ is the free space wavenumber. The geometry of the model problem is shown in Figure 1. The measurements are taken on $\mathrm{\Gamma }:=\{\mathit{x}=(x_1,x_2):x_2=b>0,-L<x_1<L\}$, and $\mathrm{\Omega }\subset {\mathbb{R}}_{+}^2$ is a compact set containing $D$ and $\mathrm{\Gamma }$.

Figure 2. Direct discretization method and Tikhonov regularization of problem (4.14.1 $$\begin{array}{c}\hfill g(t)={\int }_0^1{e}^{-ts}f(s)ds,t\in [0,1],\end{array}$$4.1 ).

Figure 3. Moment method and moment method with truncated SVD for problem (4.14.1 $$\begin{array}{c}\hfill g(t)={\int }_0^1{e}^{-ts}f(s)ds,t\in [0,1],\end{array}$$4.1 ).

Figure 4. Computing model.

Figure 5. $2c=\mathit{\lambda }/2$, $b=\mathit{\lambda }/2$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-8}$.

Figure 6. $2c=\mathit{\lambda }/2$, $b=\mathit{\lambda }/2$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-6}$.

Figure 7. $2c=\mathit{\lambda }/2$, $b=\mathit{\lambda }/2$, $\mathit{\sigma }={10}^{-6}$, $\mathrm{\Lambda }={10}^{-6}$.

Figure 8. $2c=\mathit{\lambda }/2$, $b=\mathit{\lambda }/2$, $\mathit{\sigma }={10}^{-6}$, $\mathrm{\Lambda }={10}^{-3}$.

It is easy to check that the total field $u$ satisfies the Helmholtz equation2.1 $$\begin{array}{c}\hfill \mathrm{\Delta }u+n^2{k}_0^2(1+q)u=0,\text{in}{\mathbb{R}}^2.\end{array}$$2.1 Let $u^{ref}$ be the reference field, which satisfies the Helmholtz equation without the scatterer:2.2 $$\begin{array}{c}\hfill \mathrm{\Delta }{u}^{ref}+n^2{k}_0^2{u}^{ref}=0,\text{in}{\mathbb{R}}^2.\end{array}$$2.2 From Equation (2.22.2 $$\begin{array}{c}\hfill \mathrm{\Delta }{u}^{ref}+n^2{k}_0^2{u}^{ref}=0,\text{in}{\mathbb{R}}^2.\end{array}$$2.2 ) and the continuity condition on the interface we can analytically obtain that$$u^{ref}=\left\{\begin{array}{cc}{u}^t,\hfill & \text{for}{x}_2>0,\hfill \\ u^i+u^r,\hfill & \text{for}{x}_2<0,\hfill \end{array}\right.$$where$$u^t=C_T{e}^{\mathrm{i}\mathit{\alpha }{x}_1+\mathrm{i}\mathit{\gamma }{x}_2}\text{and}{u}^r=C_R{e}^{\mathrm{i}\mathit{\alpha }{x}_1-\mathrm{i}\mathit{\beta }{x}_2}$$are the transmitted and reflected waves, respectively, [8, 12, 13], and$$C_T=\frac{2\mathit{\beta }}{\mathit{\beta }+\mathit{\gamma }},C_R=\frac{\mathit{\beta }-\mathit{\gamma }}{\mathit{\beta }+\mathit{\gamma }},$$$$\mathit{\gamma }(\mathit{\alpha })=\left\{\begin{array}{cc}\sqrt{k_0^2-{\mathit{\alpha }}^2},\hfill & \text{for}{k}_0>|\mathit{\alpha }|,\hfill \\ \mathrm{i}\sqrt{{\mathit{\alpha }}^2-k_0^2},\hfill & \text{for}{k}_0<|\mathit{\alpha }|.\hfill \end{array}\right.$$It can be readily seen that the transmitted wave $u^t$ becomes evanescent wave when $|\mathit{\alpha }|>k_0$. That is, $u^t$ decays exponentially in the $x_2$ direction.

Figure 9. $2c=\mathit{\lambda }/2$, $b=\mathit{\lambda }/2$, $\mathit{\sigma }={10}^{-3}$, $\mathrm{\Lambda }={10}^{-6}$.

Figure 10. $2c=\mathit{\lambda }/2$, $b=\mathit{\lambda }/2$, $\mathit{\sigma }={10}^{-3}$, $\mathrm{\Lambda }={10}^{-3}$.

Define the scattering field $u^s=u-u^{ref}$. Then from (2.12.1 $$\begin{array}{c}\hfill \mathrm{\Delta }u+n^2{k}_0^2(1+q)u=0,\text{in}{\mathbb{R}}^2.\end{array}$$2.1 ) and (2.22.2 $$\begin{array}{c}\hfill \mathrm{\Delta }{u}^{ref}+n^2{k}_0^2{u}^{ref}=0,\text{in}{\mathbb{R}}^2.\end{array}$$2.2 ), we have2.3 $$\begin{array}{c}\hfill \mathrm{\Delta }{u}^s+n^2{k}_0^2(1+q)u^s=-k_0^{2}q{u}^t,\text{in}{\mathbb{R}}_{+}^2.\end{array}$$2.3 And we require the radiation condition [8, 14]2.4 $$\begin{array}{c}\hfill \underset{R\to \infty }{lim}{\int }_{{\mathrm{\Sigma }}_R}\left|\frac{\mathit{\partial }{u}^s}{\mathit{\partial }\mathit{\nu }},-,\mathrm{i},n,k_0,u^s\right|dS=0,\end{array}$$2.4 where ${\mathrm{\Sigma }}_R$ is the sphere centred in the origin with radius $R$, and $\mathit{\nu }$ is the unit outward normal to ${\mathrm{\Sigma }}_R$.

Now we introduce the fundamental solution $G(\mathit{x},\mathit{y})$ of the Helmholtz equation in a two-layered medium in ${\mathbb{R}}^2$, which satisfies2.5 $$\begin{array}{c}\hfill \mathrm{\Delta }G(\mathit{x},\mathit{y})+n^2(x)k_0^{2}G(\mathit{x},\mathit{y})=-\mathit{\delta }(\mathit{x}-\mathit{y}),\end{array}$$2.5 with radiation condition, and the continuity conditions2.6 $$\begin{array}{cc}& {G(\mathit{x},\mathit{y})|}_{x_2=0^{+}}{=G(\mathit{x},\mathit{y})|}_{x_2=0^{-}},\hfill \end{array}$$2.6 2.7 $$\begin{array}{cc}& {\left.{\displaystyle \frac{\mathit{\partial }G(\mathit{x},\mathit{y})}{\mathit{\partial }{x}_2}}\right|}_{x_2=0^{+}}={\left.{\displaystyle \frac{\mathit{\partial }G(\mathit{x},\mathit{y})}{\mathit{\partial }{x}_2}}\right|}_{x_2=0^{-}}.\hfill \end{array}$$2.7 By use of fundamental solution, we can give the formula2.8 $$\begin{array}{c}\hfill u^s(\mathit{x})=k_0^2{\int }_{D}G(\mathit{x},\mathit{y})q(\mathit{y})(u^t(\mathit{y})+u^s(\mathit{y}))d\mathit{y},\end{array}$$2.8 which is well known as Lipmann-Schwinger equation.

Furthermore, if we define ${\mathit{\beta }}_1(\mathit{\xi })$ and ${\mathit{\beta }}_2(\mathit{\xi })$ by2.9 $$\begin{array}{c}\hfill {\mathit{\beta }}_1(\mathit{\xi })=\left\{\begin{array}{cc}\sqrt{k_0^2-{\mathit{\xi }}^2},\hfill & \text{for}|k_0|>\mathit{\xi },\hfill \\ \mathrm{i}\sqrt{{\mathit{\xi }}^2-k_0^2},\hfill & \text{for}|k_0|<\mathit{\xi },\hfill \end{array}\right.\end{array}$$2.9 and2.10 $$\begin{array}{c}\hfill {\mathit{\beta }}_2(\mathit{\xi })=\left\{\begin{array}{cc}\sqrt{{(n_0{k}_0)}^2-{\mathit{\xi }}^2},\hfill & \text{for}|n_0{k}_0|>\mathit{\xi },\hfill \\ \mathrm{i}\sqrt{{\mathit{\xi }}^2-{(n_0{k}_0)}^2},\hfill & \text{for}|n_0{k}_0|<\mathit{\xi },\hfill \end{array}\right.\end{array}$$2.10 respectively, then we have that2.11 $$\begin{array}{c}\hfill G(\mathit{x},\mathit{y})=\frac{\mathrm{i}}{4\mathit{\pi }}{\int }_{-\infty }^{\infty }\frac{1}{{\mathit{\beta }}_1}\left[\frac{{\mathit{\beta }}_1-{\mathit{\beta }}_2}{{\mathit{\beta }}_1+{\mathit{\beta }}_2},e^{\mathrm{i}{\mathit{\beta }}_1(x_2+y_2)},+,e^{\mathrm{i}{\mathit{\beta }}_1|x_2-y_2|}\right]e^{\mathrm{i}\mathit{\xi }(x_1-y_1)}d\mathit{\xi },\end{array}$$2.11 when $x_2>0$ and $y_2>0$.(see [8])

We rewrite (2.82.8 $$\begin{array}{c}\hfill u^s(\mathit{x})=k_0^2{\int }_{D}G(\mathit{x},\mathit{y})q(\mathit{y})(u^t(\mathit{y})+u^s(\mathit{y}))d\mathit{y},\end{array}$$2.8 ) in an operator form:2.12 $$\begin{array}{c}\hfill \left(I,-,k_0^2,T\right)u^s=k_0^{2}T{u}^t,\end{array}$$2.12 where the operator $T$ is defined as follows:2.13 $$\begin{array}{c}\hfill Tf(\mathit{x})={\int }_{D}G(\mathit{x},\mathit{y})q(\mathit{y})f(\mathit{y})d\mathit{y}.\end{array}$$2.13 If $\Vert k_0^2{T\Vert }_2<1$, ${(I-k_0^{2}T)}^{-1}$ exists and is bounded in $L^2$ norm. We can express $u^s$ by the Neumann series2.14 $$\begin{array}{c}\hfill u^s=\sum _{j=1}^{\infty }{(k_0^{2}T)}^j{u}^t.\end{array}$$2.14

Remark 1

In fact, $G(\mathit{x},\mathit{y})$ here can be split into two parts:$$G(\mathit{x},\mathit{y})=\mathrm{\Psi }(\mathit{x},\mathit{y})+\frac{\mathrm{i}}{4}{\mathrm{H}}_0^{(1)}(k_0|\mathit{x}-\mathit{y}|),$$where$$\mathrm{\Psi }(\mathit{x},\mathit{y})=\frac{\mathrm{i}}{4\mathit{\pi }}{\int }_{-\infty }^{\infty }\frac{1}{{\mathit{\beta }}_1}\left[\frac{{\mathit{\beta }}_1-{\mathit{\beta }}_2}{{\mathit{\beta }}_1+{\mathit{\beta }}_2},e^{\mathrm{i}{\mathit{\beta }}_1(x_2+y_2)}\right]e^{\mathrm{i}\mathit{\xi }(x_1-y_1)}d\mathit{\xi }.$$It is easy to see that $\mathrm{\Psi }(\mathit{x},\mathit{y})\in C^{\infty }({\mathbb{R}}_{+}^2\times {\mathbb{R}}_{+}^2)$.[15]

In [16], the authors proved that, for$$K_{\mathrm{\Omega }}f(\mathit{x}):={\int }_{\mathrm{\Omega }}\frac{\mathrm{i}}{4}{\mathrm{H}}_0^{(1)}(k_0|\mathit{x}-\mathit{y}|)f(\mathit{y})d\mathit{y},\forall \mathit{x}\in {\mathbb{R}}^2,$$$K_{\mathrm{\Omega }}$ is a bounded operator from $L^2(\mathrm{\Omega })$ to $L^2(\mathrm{\Omega })$.

So, the operator $T$ is also bounded from $L^2(\mathrm{\Omega })$ to $L^2(\mathrm{\Omega })$. And when $k_0$ and ${\Vert q\Vert }_{L^2(D)}$ are sufficiently small, the assumption $\Vert k_0^{2}T\Vert ⩽1$ is valid.

2.2 Inverse scattering problem

The inverse scattering problem considered in this paper is from some model of near-field optics. In these problems, the measurements of scattering wave $u^s$ are given on $\mathrm{\Gamma }:=\{\mathit{x}=(x_1,x_2):x_2=b>0,-L<x_1<L\}$. The inverse problem is to determine $q$. In particular, the shape of support of $q$ is easy to be seen after that $q$ is determined.

Because the operator $T$ is linearly in $q$, our inverse problem can be written as: finding $q$ to satisfy2.15 $$\begin{array}{c}\hfill \mathit{\varphi }=K_{1}q+K_{2}q\otimes q+K_{3}q\otimes q\otimes q+\cdots ,\end{array}$$2.15 where $\mathit{\varphi }=u^s{|}_{\mathrm{\Gamma }}$ and $K_j:L^2(D)\times \cdots \times L^2(D)\to L^2(\mathrm{\Gamma })$ defined by$$\begin{array}{cc}\hfill K_{1}q& =[k_0^{2}T{u}^t]{|}_{\mathrm{\Gamma }},\hfill \\ \hfill K_{2}q\otimes q& =[{(k_0^{2}T)}^2{u}^t]{|}_{\mathrm{\Gamma }},\hfill \\ & \cdots \hfill \end{array}$$and so on.

For $q\in L^2(D)$, we have that2.16 $$\begin{array}{c}\hfill \Vert K_j{\Vert }_2⩽{\mathit{\nu }}_2{\mathit{\mu }}_2^{j-1},\end{array}$$2.16 where2.17 $$\begin{array}{c}\hfill {\mathit{\mu }}_2=k_0^2\underset{\mathit{x}\in \mathrm{\Gamma }}{sup}{\Vert G(\mathit{x},·)\Vert }_{L^2(D)}\end{array}$$2.17 and2.18 $$\begin{array}{c}\hfill {\mathit{\nu }}_2=k_0^2\Vert u^t{\Vert }_{L^{\infty }(D)}\underset{\mathit{x}\in \mathrm{\Gamma }}{sup}{\Vert G(\mathit{x},·)\Vert }_{L^2(D)}.\end{array}$$2.18 By Denoting $K_1^{+}$ as the regularized pseudo-inverse of the operator $K_1$, we have the following result from [11] and [17]:

Theorem 2.1

Suppose that $\Vert K_1^{+}{\Vert }_2<1/({\mathit{\mu }}_2+{\mathit{\nu }}_2),{\Vert K_1^{+}\mathit{\varphi }\Vert }_{L^2(D)}<1/({\mathit{\mu }}_2+{\mathit{\nu }}_2)$. Let $\mathcal{M}=max\{\Vert q{\Vert }_{L^p(D)},\Vert K_1^{+}{K}_{1}q{\Vert }_{L^p(D)}\}$ and assume that $\mathcal{M}<1/({\mathit{\mu }}_2+{\mathit{\nu }}_2)$. Then the norm of the difference between the true scatterer and the linearized inverse solution obeys the estimate2.19 $$\begin{array}{c}\hfill {∥q,-,K_1^{+},\mathit{\varphi }∥}_{L^2(D)}⩽C\frac{[({\mathit{\mu }}_2+{\mathit{\nu }}_2)\Vert K_1^{+}\mathit{\varphi }{\Vert }_2]{}^2}{1-({\mathit{\mu }}_2+{\mathit{\nu }}_2){\Vert K_1^{+}\mathit{\varphi }\Vert }_2}+\tilde{C}{\Vert (I-K_1^{+}{K}_1)q\Vert }_{L^2(D)},\end{array}$$2.19 where $C$ and $\tilde{C}$ are independent of $\mathit{\varphi }$.

Theorem 2.1 gives the error estimate of the inverse problem concerning Born approximation. The first error term in (2.192.19 $$\begin{array}{c}\hfill {∥q,-,K_1^{+},\mathit{\varphi }∥}_{L^2(D)}⩽C\frac{[({\mathit{\mu }}_2+{\mathit{\nu }}_2)\Vert K_1^{+}\mathit{\varphi }{\Vert }_2]{}^2}{1-({\mathit{\mu }}_2+{\mathit{\nu }}_2){\Vert K_1^{+}\mathit{\varphi }\Vert }_2}+\tilde{C}{\Vert (I-K_1^{+}{K}_1)q\Vert }_{L^2(D)},\end{array}$$2.19 ) relates much to the the norm of $K_1^{+}\mathit{\varphi }$. And the second term describes the error between the exact solution and the regularized solution of the linearized problem with $K_{1}q$ being the measured data.

So the construction of $K_1^{+}$ is very important to the inverse scattering problem. However, the equation $K_{1}q=\mathit{\varphi }$ is difficult to solve, for its severe ill-posedness. Few work has been done for the detailed numerical methods of this severely ill-posed problem and its error estimate, which we will discuss in the next section.

Figure 11. $2c=\mathit{\lambda }/2$, $b=\mathit{\lambda }$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-6}$.

Figure 12. $2c=\mathit{\lambda }/2$, $b=\mathit{\lambda }$, $\mathit{\sigma }={10}^{-6}$, $\mathrm{\Lambda }={10}^{-6}$.

Figure 13. $2c=\mathit{\lambda }/2$, $b=\mathit{\lambda }$, $\mathit{\sigma }={10}^{-3}$, $\mathrm{\Lambda }={10}^{-3}$.

Figure 14. $2c=\mathit{\lambda }/4$, $b=\mathit{\lambda }/2$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-6}$.

Remark 2

When $k_0$, $\mathrm{s}upp(q)$ and ${\Vert q\Vert }_{L^2(D)}$ are sufficiently small, the assumption $\Vert k_0^{2}T\Vert ⩽1$ is valid, and ${\mathit{\mu }}_2,$${\mathit{\nu }}_2$ are also sufficiently small.

We checked numerically that $({\mathit{\mu }}_2+{\mathit{\nu }}_2){\Vert K_1^{+}\Vert }_2$, $({\mathit{\mu }}_2+{\mathit{\nu }}_2){\Vert K_1^{+}\mathit{\varphi }\Vert }_2$ are very small. Therefore, the assumptions in Theorem 2.1 could be satisfied, and the first term in (2.192.19 $$\begin{array}{c}\hfill {∥q,-,K_1^{+},\mathit{\varphi }∥}_{L^2(D)}⩽C\frac{[({\mathit{\mu }}_2+{\mathit{\nu }}_2)\Vert K_1^{+}\mathit{\varphi }{\Vert }_2]{}^2}{1-({\mathit{\mu }}_2+{\mathit{\nu }}_2){\Vert K_1^{+}\mathit{\varphi }\Vert }_2}+\tilde{C}{\Vert (I-K_1^{+}{K}_1)q\Vert }_{L^2(D)},\end{array}$$2.19 ) would also be small.

3 Numerical method for the linearized inverse problem

The linearized inverse problem, as discussed in pervious section, is3.1 $$\begin{array}{c}\hfill u^s=k_0^{2}T{u}^t,\end{array}$$3.1 more explicitly,3.2 $$\begin{array}{c}\hfill u^s(\mathit{x})=k_0^2{\int }_{D}G(\mathit{x},\mathit{y})q(\mathit{y})u^t(\mathit{y})d\mathit{y}.\end{array}$$3.2 Note that the only known data lie on the segment $\mathrm{\Gamma }=\{\mathit{x}:x_2=b,-L<x_1<L\}$, so we extend $u^s$ on $\mathrm{\Gamma }$ to the whole line $\{\mathit{x}:x_2=b,-\infty <x_1<\infty \}$ by defining $u^s=0$ on the remaining interval of $\mathrm{\Gamma }$. We assume that the error between the expanded data $u^s$ and the “exact” data is small when $L$ is sufficiently large. Namely, we assume that $\Vert u^s-u_{\mathrm{e}xact}^s\Vert <{\mathit{\delta }}_1$, where ${\mathit{\delta }}_1$ depends on $L$. The scattering data $u^s$ depends on $\mathit{\alpha }$, and we need several different incident waves, i.e. we need $[{\mathit{\alpha }}_1,\cdots ,{\mathit{\alpha }}_n]\in [-n_0{k}_0,n_0{k}_0]$ to reconstruct the scatterer $q(\mathit{x})$.

Thus, we have3.3 $$\begin{array}{c}\hfill u^s(x_1,b)=k_0^2{\int }_{-\infty }^{\infty }{\int }_0^{b}G(x_1,b;y_1,y_2)q(y_1,y_2)u^t(y_1,y_2)d{y}_{2}d{y}_1.\end{array}$$3.3 Substituting the fundamental solution $G(\mathit{x},\mathit{y})$ and $u^t$ into (3.33.3 $$\begin{array}{c}\hfill u^s(x_1,b)=k_0^2{\int }_{-\infty }^{\infty }{\int }_0^{b}G(x_1,b;y_1,y_2)q(y_1,y_2)u^t(y_1,y_2)d{y}_{2}d{y}_1.\end{array}$$3.3 ) and realizing the dependence of $u^s$ on $\mathit{\alpha }$, we rewrite (3.33.3 $$\begin{array}{c}\hfill u^s(x_1,b)=k_0^2{\int }_{-\infty }^{\infty }{\int }_0^{b}G(x_1,b;y_1,y_2)q(y_1,y_2)u^t(y_1,y_2)d{y}_{2}d{y}_1.\end{array}$$3.3 ) as3.4 $$\begin{array}{cc}\hfill u^s(x_1,b;\mathit{\alpha })& =\frac{\mathrm{i}{k}_0^2{C}_T}{4\mathit{\pi }}{\int }_{-\infty }^{\infty }{\int }_0^{b}q(y_1,y_2)e^{\mathrm{i}\mathit{\alpha }{y}_1+\mathrm{i}\mathit{\gamma }{y}_2}\hfill \\ & \times {\int }_{-\infty }^{\infty }\frac{1}{{\mathit{\beta }}_1}\left[\frac{{\mathit{\beta }}_1-{\mathit{\beta }}_2}{{\mathit{\beta }}_1+{\mathit{\beta }}_2},e^{\mathrm{i}{\mathit{\beta }}_1(b+y_2)},+,e^{\mathrm{i}{\mathit{\beta }}_1(b-y_2)}\right]e^{\mathrm{i}\mathit{\xi }(x_1-y_1)}d\mathit{\xi }d{y}_{2}d{y}_1.\hfill \end{array}$$3.4 Multiplying both sides of (3.43.4 $$\begin{array}{cc}\hfill u^s(x_1,b;\mathit{\alpha })& =\frac{\mathrm{i}{k}_0^2{C}_T}{4\mathit{\pi }}{\int }_{-\infty }^{\infty }{\int }_0^{b}q(y_1,y_2)e^{\mathrm{i}\mathit{\alpha }{y}_1+\mathrm{i}\mathit{\gamma }{y}_2}\hfill \\ & \times {\int }_{-\infty }^{\infty }\frac{1}{{\mathit{\beta }}_1}\left[\frac{{\mathit{\beta }}_1-{\mathit{\beta }}_2}{{\mathit{\beta }}_1+{\mathit{\beta }}_2},e^{\mathrm{i}{\mathit{\beta }}_1(b+y_2)},+,e^{\mathrm{i}{\mathit{\beta }}_1(b-y_2)}\right]e^{\mathrm{i}\mathit{\xi }(x_1-y_1)}d\mathit{\xi }d{y}_{2}d{y}_1.\hfill \end{array}$$3.4 ) by $e^{-\mathrm{i}\mathit{\omega }{x}_1}$ and integrating from $-\infty$ to $\infty$ with respect to $x_1$,3.5 $$\begin{array}{cc}& {\widehat{u}}^s(\mathit{\omega },b;\mathit{\alpha })=\frac{\mathrm{i}{k}_0^2{C}_T}{4\mathit{\pi }}{\int }_{-\infty }^{\infty }{\int }_0^{b}q(y_1,y_2)e^{\mathrm{i}\mathit{\alpha }{y}_1+\mathrm{i}\mathit{\gamma }{y}_2}\hfill \\ & \times {\int }_{-\infty }^{\infty }\frac{1}{{\mathit{\beta }}_1}\left[\frac{{\mathit{\beta }}_1-{\mathit{\beta }}_2}{{\mathit{\beta }}_1+{\mathit{\beta }}_2},e^{\mathrm{i}{\mathit{\beta }}_1(b+y_2)},+,e^{\mathrm{i}{\mathit{\beta }}_1(b-y_2)}\right]e^{-\mathrm{i}\mathit{\xi }{y}_1}{\int }_{-\infty }^{\infty }{e}^{\mathrm{i}(\mathit{\xi }-\mathit{\omega })x_1}d{x}_{1}d\mathit{\xi }d{y}_{2}d{y}_1,\hfill \end{array}$$3.5 and noting that$${\int }_{-\infty }^{\infty }{e}^{\mathrm{i}(\mathit{\xi }-\mathit{\omega })x_1}d{x}_1=2\mathit{\pi }\mathit{\delta }(\mathit{\xi }-\mathit{\omega }),$$(3.53.5 $$\begin{array}{cc}& {\widehat{u}}^s(\mathit{\omega },b;\mathit{\alpha })=\frac{\mathrm{i}{k}_0^2{C}_T}{4\mathit{\pi }}{\int }_{-\infty }^{\infty }{\int }_0^{b}q(y_1,y_2)e^{\mathrm{i}\mathit{\alpha }{y}_1+\mathrm{i}\mathit{\gamma }{y}_2}\hfill \\ & \times {\int }_{-\infty }^{\infty }\frac{1}{{\mathit{\beta }}_1}\left[\frac{{\mathit{\beta }}_1-{\mathit{\beta }}_2}{{\mathit{\beta }}_1+{\mathit{\beta }}_2},e^{\mathrm{i}{\mathit{\beta }}_1(b+y_2)},+,e^{\mathrm{i}{\mathit{\beta }}_1(b-y_2)}\right]e^{-\mathrm{i}\mathit{\xi }{y}_1}{\int }_{-\infty }^{\infty }{e}^{\mathrm{i}(\mathit{\xi }-\mathit{\omega })x_1}d{x}_{1}d\mathit{\xi }d{y}_{2}d{y}_1,\hfill \end{array}$$3.5 ) can be rewritten as3.6 $$\begin{array}{cc}\hfill {\widehat{u}}^s(\mathit{\omega },b;\mathit{\alpha })& =\frac{\mathrm{i}{k}_0^2{C}_T}{2}{\int }_{-\infty }^{\infty }{\int }_0^{b}q(y_1,y_2)e^{\mathrm{i}\mathit{\alpha }{y}_1+\mathrm{i}\mathit{\gamma }{y}_2}\hfill \\ & ·\frac{1}{{\mathit{\beta }}_1(\mathit{\omega })}\left[\frac{{\mathit{\beta }}_1(\mathit{\omega })-{\mathit{\beta }}_2(\mathit{\omega })}{{\mathit{\beta }}_1(\mathit{\omega })+{\mathit{\beta }}_2(\mathit{\omega })},e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })(b+y_2)},+,e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })(b-y_2)}\right]e^{-\mathrm{i}\mathit{\omega }{y}_1}d{y}_{2}d{y}_1.\hfill \end{array}$$3.6 Realizing that$${\int }_{-\infty }^{\infty }q(y_1,y_2)e^{\mathrm{i}\mathit{\alpha }{y}_1}{e}^{-\mathrm{i}\mathit{\omega }{y}_1}d{y}_1=\widehat{q}(\mathit{\omega }-\mathit{\alpha },y_2),$$(3.63.6 $$\begin{array}{cc}\hfill {\widehat{u}}^s(\mathit{\omega },b;\mathit{\alpha })& =\frac{\mathrm{i}{k}_0^2{C}_T}{2}{\int }_{-\infty }^{\infty }{\int }_0^{b}q(y_1,y_2)e^{\mathrm{i}\mathit{\alpha }{y}_1+\mathrm{i}\mathit{\gamma }{y}_2}\hfill \\ & ·\frac{1}{{\mathit{\beta }}_1(\mathit{\omega })}\left[\frac{{\mathit{\beta }}_1(\mathit{\omega })-{\mathit{\beta }}_2(\mathit{\omega })}{{\mathit{\beta }}_1(\mathit{\omega })+{\mathit{\beta }}_2(\mathit{\omega })},e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })(b+y_2)},+,e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })(b-y_2)}\right]e^{-\mathrm{i}\mathit{\omega }{y}_1}d{y}_{2}d{y}_1.\hfill \end{array}$$3.6 ) can be reduced to3.7 $$\begin{array}{cc}& {\widehat{u}}^s(\mathit{\omega },b;\mathit{\alpha })\hfill \\ & =\frac{\mathrm{i}{k}_0^2{C}_T}{2}{\int }_0^b\frac{e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })b}}{{\mathit{\beta }}_1(\mathit{\omega })}\left[\frac{{\mathit{\beta }}_1(\mathit{\omega })-{\mathit{\beta }}_2(\mathit{\omega })}{{\mathit{\beta }}_1(\mathit{\omega })+{\mathit{\beta }}_2(\mathit{\omega })},e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2},+,e^{-\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2}\right]e^{\mathrm{i}\mathit{\gamma }{y}_2}\widehat{q}(\mathit{\omega }-\mathit{\alpha },y_2)d{y}_2.\hfill \end{array}$$3.7 When $|\mathit{\omega }|<k_0$ and $|\mathit{\alpha }|<k_0$, both ${\mathit{\beta }}_1$ and $\mathit{\gamma }$ are real; and for $k_0<|\mathit{\alpha }|<n_{0}k$, $\mathit{\gamma }$ is pure imaginary, and so is ${\mathit{\beta }}_1$ with $|\mathit{\omega }|>k_0$. It is obvious that the inversion of (3.73.7 $$\begin{array}{cc}& {\widehat{u}}^s(\mathit{\omega },b;\mathit{\alpha })\hfill \\ & =\frac{\mathrm{i}{k}_0^2{C}_T}{2}{\int }_0^b\frac{e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })b}}{{\mathit{\beta }}_1(\mathit{\omega })}\left[\frac{{\mathit{\beta }}_1(\mathit{\omega })-{\mathit{\beta }}_2(\mathit{\omega })}{{\mathit{\beta }}_1(\mathit{\omega })+{\mathit{\beta }}_2(\mathit{\omega })},e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2},+,e^{-\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2}\right]e^{\mathrm{i}\mathit{\gamma }{y}_2}\widehat{q}(\mathit{\omega }-\mathit{\alpha },y_2)d{y}_2.\hfill \end{array}$$3.7 ) involves both inverse Fourier transform and inverse Laplace transform. In fact, both inverse Fourier transform and inverse Laplace transform defined on finite interval are ill-posed, especially the inverse Laplace transform is exponentially ill-posed.[18]

In [9], the authors used a method that is essentially the moment method to handle the inversion of an optical diffusion problem. They simulated their problem in the two-dimensional case, but the resolution is limited by numerical precisions.

Here, we split the two-dimensional problem (3.73.7 $$\begin{array}{cc}& {\widehat{u}}^s(\mathit{\omega },b;\mathit{\alpha })\hfill \\ & =\frac{\mathrm{i}{k}_0^2{C}_T}{2}{\int }_0^b\frac{e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })b}}{{\mathit{\beta }}_1(\mathit{\omega })}\left[\frac{{\mathit{\beta }}_1(\mathit{\omega })-{\mathit{\beta }}_2(\mathit{\omega })}{{\mathit{\beta }}_1(\mathit{\omega })+{\mathit{\beta }}_2(\mathit{\omega })},e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2},+,e^{-\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2}\right]e^{\mathrm{i}\mathit{\gamma }{y}_2}\widehat{q}(\mathit{\omega }-\mathit{\alpha },y_2)d{y}_2.\hfill \end{array}$$3.7 ) into two one-dimensional problems. So, the computation runs much faster and uses less memory. Moreover, we apply the moment method with truncated SVD to get a better resolution where the limitation of numerical precisions is circumvented in a certain extent.

Figure 15. $2c=\mathit{\lambda }/4$, $b=\mathit{\lambda }/2$, $\mathit{\sigma }={10}^{-3}$, $\mathrm{\Lambda }={10}^{-3}$.

We need two steps to solve (3.73.7 $$\begin{array}{cc}& {\widehat{u}}^s(\mathit{\omega },b;\mathit{\alpha })\hfill \\ & =\frac{\mathrm{i}{k}_0^2{C}_T}{2}{\int }_0^b\frac{e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })b}}{{\mathit{\beta }}_1(\mathit{\omega })}\left[\frac{{\mathit{\beta }}_1(\mathit{\omega })-{\mathit{\beta }}_2(\mathit{\omega })}{{\mathit{\beta }}_1(\mathit{\omega })+{\mathit{\beta }}_2(\mathit{\omega })},e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2},+,e^{-\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2}\right]e^{\mathrm{i}\mathit{\gamma }{y}_2}\widehat{q}(\mathit{\omega }-\mathit{\alpha },y_2)d{y}_2.\hfill \end{array}$$3.7 ).

Step 1

For $\mathit{\omega }-\mathit{\alpha }=c^{(k)},k=1,2,\cdots ,$ obtain $\widehat{q}(c^{(k)},y_2)$ from a one-dimensional integral equation for every $k$;

So for every calculation of $\widehat{q}(c^{(k)},y_2),k=1,2,\cdots$, we need special choices of $\mathit{\omega }$. That is, if we have $n$ incident angles, then we set ${\mathit{\omega }}_j={\mathit{\alpha }}_j+c^{(k)},j=1,\cdots ,n,$ for $k=1,2,\cdots$. We suppose all the ${\mathit{\omega }}_j$’s lie in $[-W,W]$.

Step 2

Obtain $q(y_1,y_2)$ from $\widehat{q}(c^{(k)},y_2)$ by inverse Fourier transform on finite intervals.

After the two steps, an approximate solution of the linearized inverse problem (3.13.1 $$\begin{array}{c}\hfill u^s=k_0^{2}T{u}^t,\end{array}$$3.1 ) is obtained. Obviously, the computation of step 1 is quite difficult for its severe ill-posedness, especially on a finite data-set. As is discussed above, we treat it as a kind of perturbation of the measured data.

Figure 16. $2c=\mathit{\lambda }/4$, $b=\mathit{\lambda }$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-6}$.

Figure 17. $2c=\mathit{\lambda }/4$, $b=\mathit{\lambda }$, $\mathit{\sigma }={10}^{-3}$, $\mathrm{\Lambda }={10}^{-3}$.

The first step is to solve the operator equation for every specified $k,k=1,2,\cdots$.3.8 $$\begin{array}{c}\hfill (K\widehat{q})(\mathit{\omega })={\widehat{u}}^s(\mathit{\omega },\mathit{\alpha })\frac{2}{\mathrm{i}{k}_0^2{C}_T}{\mathit{\beta }}_1(\mathit{\omega }),\end{array}$$3.8 where3.9 $$\begin{array}{cc}& (K\widehat{q})(\mathit{\omega })\hfill \\ & :={\int }_0^b{e}^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })b}\left[\frac{{\mathit{\beta }}_1(\mathit{\omega })-{\mathit{\beta }}_2(\mathit{\omega })}{{\mathit{\beta }}_1(\mathit{\omega })+{\mathit{\beta }}_2(\mathit{\omega })},e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2},+,e^{-\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2}\right]e^{\mathrm{i}\mathit{\gamma }(\mathit{\omega }-c^{(k)})y_2}\widehat{q}(c^{(k)},y_2)d{y}_2.\hfill \end{array}$$3.9 Denote$$k(\mathit{\omega },y_2)=e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })b}\left[\frac{{\mathit{\beta }}_1(\mathit{\omega })-{\mathit{\beta }}_2(\mathit{\omega })}{{\mathit{\beta }}_1(\mathit{\omega })+{\mathit{\beta }}_2(\mathit{\omega })},e^{\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2},+,e^{-\mathrm{i}{\mathit{\beta }}_1(\mathit{\omega })y_2}\right]e^{\mathrm{i}\mathit{\gamma }(\mathit{\omega }-c^{(k)})y_2},$$then $K:L^2[0,b]\to L^2[-W,W]$ is a compact operator since $k(\mathit{\omega },y_2)\in L^2([-W,W]\times [0,b])$.

Let $X_n\subset L^2[0,b]$ be finite-dimensional subspaces with $dim{X}_n=n,$ and $0⩽{\mathit{\omega }}_1<\cdots <{\mathit{\omega }}_n⩽b$ be the collocation points, which are chosen as ${\mathit{\omega }}_j={\mathit{\alpha }}_j+c^{(k)},j=1,\cdots ,n$. Then the collocation equation is3.10 $$\begin{array}{c}\hfill (K{\widehat{q}}_n)({\mathit{\omega }}_i)=y({\mathit{\omega }}_i),i=1,...,n,\end{array}$$3.10 where$$y({\mathit{\omega }}_i)={\widehat{u}}^s({\mathit{\omega }}_i,\mathit{\alpha })\frac{2}{\mathrm{i}{k}_0^2{C}_T}{\mathit{\beta }}_1({\mathit{\omega }}_i),$$then ${\widehat{q}}_n\in X_n$ is a collocation solution of (3.103.10 $$\begin{array}{c}\hfill (K{\widehat{q}}_n)({\mathit{\omega }}_i)=y({\mathit{\omega }}_i),i=1,...,n,\end{array}$$3.10 ).

Choosing a basis $\{{\mathit{\varphi }}_j,j=1,\cdots ,n\}$ of $X_n,$ then ${\widehat{q}}_n(c^{(k)},y_2)$ can be represented as$${\widehat{q}}_n(c^{(k)},y_2)=\sum _{j=1}^n{a}_j{\mathit{\varphi }}_j(y_2),$$where the coefficients $a_j$’s are the solution of the system3.11 $$\begin{array}{c}\hfill Aa=\mathit{\beta }\end{array}$$3.11 with$$A_{i,j}=K{\mathit{\varphi }}_j({\mathit{\omega }}_i),{\mathit{\beta }}_i=y({\mathit{\omega }}_i).$$For ill-posed problems such as (3.113.11 $$\begin{array}{c}\hfill Aa=\mathit{\beta }\end{array}$$3.11 ), even small perturbation of data $\mathit{\beta }$ can change the solution a lot. Recall that we can only obtain the scattering data on $\mathrm{\Gamma }$, and we have assumed that $\Vert u^s-u_{\mathrm{e}xact}^s\Vert <{\mathit{\delta }}_1$. Now we denote the perturbed data, or the measured data, by ${\mathit{\beta }}^{\mathit{\delta }}$, and let $|\mathit{\beta }-{\mathit{\beta }}^{\mathit{\delta }}|<\mathit{\delta }$, where $\mathit{\delta }$ takes in the influence of the error bound ${\mathit{\delta }}_1$, and $|·|$ denotes Euclidean norm in ${\mathbb{C}}^n.$

Remark 3

In fact, there is no need that the number of measurements points coincides that of the collocation points. In [19], a proper interpolation method can be made to obtain all the $y({\mathit{\omega }}_j)$’s ,$j=1,\cdots ,n$, from the measurements points $s_j,j=1,\cdots ,m.$

Solving (3.103.10 $$\begin{array}{c}\hfill (K{\widehat{q}}_n)({\mathit{\omega }}_i)=y({\mathit{\omega }}_i),i=1,...,n,\end{array}$$3.10 ) in $L^2[0,b]$ cannot specify the solution ${\widehat{q}}_n$ uniquely. There are several choices of the basis ${\mathit{\varphi }}_j$’s. In this section, we will give an easy choice of the basis, and under this choice, we can obtain a unique solution.

We will find a moment solution of (3.103.10 $$\begin{array}{c}\hfill (K{\widehat{q}}_n)({\mathit{\omega }}_i)=y({\mathit{\omega }}_i),i=1,...,n,\end{array}$$3.10 ), which means that$$\Vert x_n{\Vert }_{L^2}=min\{\Vert z_n{\Vert }_{L^2}:z_n\in L^2[0,b],z_n\text{satisfies}(3.10)\}$$with respect to the collocation points ${\{{\mathit{\omega }}_j\}}_{j=1}^n$.

Recall Riesz theorem that there exists $k_j\in L^2[0,b]$, such that, $Kz(t_j)=(z,k_j)$, for all $z\in L^2[0,b]$, and $j=1,\cdots ,n.$ In particular for the integral operator $K$, $k_j$ is explicitly given by$$k_j(y_2)=k({\mathit{\omega }}_j,y_2),j=1,\cdots ,n.$$If $k_j$’s are linearly independent, i.e. the determinant of the collocation matrix $A_{i,j}=(k_i,k_j)$ is not equal to 0, by choosing ${\mathit{\varphi }}_j=k_j$ , there exists one and only one moment solution ${\widehat{q}}_n$ of (3.103.10 $$\begin{array}{c}\hfill (K{\widehat{q}}_n)({\mathit{\omega }}_i)=y({\mathit{\omega }}_i),i=1,...,n,\end{array}$$3.10 ) (see [20]).

However, the system of (3.113.11 $$\begin{array}{c}\hfill Aa=\mathit{\beta }\end{array}$$3.11 ) is severely ill-posed, so it is difficult to find a stable moment solution. We will apply the truncated SVD method to improve the stability of problem (3.103.10 $$\begin{array}{c}\hfill (K{\widehat{q}}_n)({\mathit{\omega }}_i)=y({\mathit{\omega }}_i),i=1,...,n,\end{array}$$3.10 ).

Using SVD, we can write$$A=U\mathrm{\Sigma }{V}^{\ast },$$where $U$ and $V$ are unitary, and$$U=[u_1,\cdots ,u_n],V=[v_1,\cdots ,v_n],\mathrm{\Sigma }=\mathrm{d}iag({\mathit{\lambda }}_1,\cdots ,{\mathit{\lambda }}_n).$$Thus,$$A{v}_j={\mathit{\lambda }}_j{u}_j,A^{\ast }{u}_j={\mathit{\lambda }}_j{v}_j.$$We assume that ${\mathit{\lambda }}_1⩾\cdots ⩾{\mathit{\lambda }}_n⩾0$, and let $\mathrm{\Lambda }⩾{\mathit{\lambda }}_K$ for some $K,1⩽K⩽n$. And denote $\tilde{\mathrm{\Sigma }}=\mathrm{d}iag({\mathit{\lambda }}_1,\cdots ,{\mathit{\lambda }}_K,0,\cdots ,0)$. Let $\tilde{A}=U\tilde{\mathrm{\Sigma }}{V}^{\ast }$, we have the following estimate.

Theorem 3.1

Let $\{{\mathit{\omega }}_1^{(n)},{\mathit{\omega }}_2^{(n)},\cdots ,{\mathit{\omega }}_n^{(n)}\}\subset [0,b],n\in \mathbb{N},$ be a sequence of collocation points and$$X_n:=\mathrm{s}pan\{k_j^{(n)},j=1,\cdots ,n\}.$$And assume that there exist $z\in {\mathbb{C}}^n$ such that $a=A^{\ast }z$. If ${\tilde{\widehat{q}}}_n^{\mathit{\delta }}={\sum }_{j=1}^n{\tilde{a}}_j^{\mathit{\delta }}{k}_j^{(n)}$, where ${\tilde{a}}^{\mathit{\delta }}\in C^n$ solves $\tilde{A}{\tilde{a}}^{\mathit{\delta }}={\mathit{\beta }}^{\mathit{\delta }}$ with $|{\mathit{\beta }}^{\mathit{\delta }}-\mathit{\beta }|⩽\mathit{\delta }$, then the following error estimate holds:3.12 $$\begin{array}{c}\hfill \Vert {\tilde{\widehat{q}}}_n^{\mathit{\delta }}-\widehat{q}\Vert ⩽p_n\left(\frac{\mathit{\delta }}{\mathrm{\Lambda }},+,\mathrm{\Lambda },|z|\right)+cmin\{\Vert \widehat{q}-z_n{\Vert }_{L^2}:z_n\in X_n\},\end{array}$$3.12 where$$p_n=max\left\{{∥\sum _{j=1}^n,{\mathit{\rho }}_j,k_j^{(n)}∥}_{L^2},:,\sum _{j=1}^n,{|{\mathit{\rho }}_j|}^2,=,1\right\}.$$

Proof

We write $\Vert {\tilde{\widehat{q}}}_n^{\mathit{\delta }}-\widehat{q}\Vert ⩽\Vert {\tilde{\widehat{q}}}_n^{\mathit{\delta }}-{\widehat{q}}_n\Vert +\Vert {\widehat{q}}_n-\widehat{q}\Vert$, and the second term is estimated in [20]. All we have to do in this proof is to estimate $\Vert {\tilde{\widehat{q}}}_n^{\mathit{\delta }}-{\widehat{q}}_n\Vert$.3.13 $$\begin{array}{cc}\hfill \Vert {\tilde{\widehat{q}}}_n^{\mathit{\delta }}-{\widehat{q}}_n\Vert & =∥\sum _{j=1}^n,{\tilde{a}}_j^{\mathit{\delta }},k_j^{(n)},-,\sum _{j=1}^n,a_j,k_j^{(n)}∥\hfill \\ & ⩽\hfill & \hfill p_n|{\tilde{a}}^{\mathit{\delta }}-a|=p_n|{\tilde{A}}^{-1}{\mathit{\beta }}^{\mathit{\delta }}-A^{-1}\mathit{\beta }|\\ & =p_n|{\tilde{A}}^{-1}{\mathit{\beta }}^{\mathit{\delta }}-{\tilde{A}}^{-1}\mathit{\beta }+{\tilde{A}}^{-1}\mathit{\beta }-A^{-1}\mathit{\beta }|\hfill \\ & ⩽\hfill & \hfill p_n(|{\tilde{A}}^{-1}({\mathit{\beta }}^{\mathit{\delta }}-\mathit{\beta })|+|({\tilde{A}}^{-1}-A^{-1})\mathit{\beta }|),\end{array}$$3.13 the first term in (3.162.16 $$\begin{array}{c}\hfill \Vert K_j{\Vert }_2⩽{\mathit{\nu }}_2{\mathit{\mu }}_2^{j-1},\end{array}$$2.16 ) is$$|{\tilde{A}}^{-1}({\mathit{\beta }}^{\mathit{\delta }}-\mathit{\beta })|⩽{\mathit{\lambda }}_K^{-1}\mathit{\delta }⩽\frac{\mathit{\delta }}{\mathrm{\Lambda }},$$and the second term in (3.162.16 $$\begin{array}{c}\hfill \Vert K_j{\Vert }_2⩽{\mathit{\nu }}_2{\mathit{\mu }}_2^{j-1},\end{array}$$2.16 ) is3.14 $$\begin{array}{cc}\hfill |({\tilde{A}}^{-1}-A^{-1})\mathit{\beta }|& =\left|\sum _{j=K+1}^n,\frac{1}{{\mathit{\lambda }}_j},(\mathit{\beta },u_j),v_j\right|\hfill \\ & =\left|\sum _{j=K+1}^n,\frac{1}{{\mathit{\lambda }}_j},(Aa,u_j),v_j\right|\hfill \\ & =\left|\sum _{j=K+1}^n,(a,v_j),v_j\right|.\hfill \end{array}$$3.14 If $a=A^{\ast }z,$ then$$(a,v_j)=(A^{\ast }z,v_j)=(z,A{v}_j)=(z,{\mathit{\lambda }}_j{u}_j)={\mathit{\lambda }}_j(a,u_j).$$Thus, (3.192.19 $$\begin{array}{c}\hfill {∥q,-,K_1^{+},\mathit{\varphi }∥}_{L^2(D)}⩽C\frac{[({\mathit{\mu }}_2+{\mathit{\nu }}_2)\Vert K_1^{+}\mathit{\varphi }{\Vert }_2]{}^2}{1-({\mathit{\mu }}_2+{\mathit{\nu }}_2){\Vert K_1^{+}\mathit{\varphi }\Vert }_2}+\tilde{C}{\Vert (I-K_1^{+}{K}_1)q\Vert }_{L^2(D)},\end{array}$$2.19 ) can be written as$$\left|\sum _{j=K+1}^n,{\mathit{\lambda }}_j,(a,u_j),v_j\right|⩽\mathrm{\Lambda }|z|.$$And the proof ends. $\square$

Without the truncation, (3.123.12 $$\begin{array}{c}\hfill \Vert {\tilde{\widehat{q}}}_n^{\mathit{\delta }}-\widehat{q}\Vert ⩽p_n\left(\frac{\mathit{\delta }}{\mathrm{\Lambda }},+,\mathrm{\Lambda },|z|\right)+cmin\{\Vert \widehat{q}-z_n{\Vert }_{L^2}:z_n\in X_n\},\end{array}$$3.12 ) is$$\Vert {\widehat{q}}_n^{\mathit{\delta }}-\widehat{q}\Vert ⩽\frac{p_n\mathit{\delta }}{{\mathit{\lambda }}_n}+cmin\{\Vert \widehat{q}-z_n{\Vert }_{L^2}:z_n\in X_n\},$$as in [20], and it is the error bound of the moment method.

Generally, we need the dimension of $X_n$ large enough to make the second term in (3.123.12 $$\begin{array}{c}\hfill \Vert {\tilde{\widehat{q}}}_n^{\mathit{\delta }}-\widehat{q}\Vert ⩽p_n\left(\frac{\mathit{\delta }}{\mathrm{\Lambda }},+,\mathrm{\Lambda },|z|\right)+cmin\{\Vert \widehat{q}-z_n{\Vert }_{L^2}:z_n\in X_n\},\end{array}$$3.12 ) to be small, but if ${\mathit{\lambda }}_j$’s decay rapidly, especially for exponentially ill-posed problems as in this paper, ${\mathit{\lambda }}_j\sim e^{-j}$, and $\frac{\mathit{\delta }}{{\mathit{\lambda }}_n}\to \infty$ as $n\to \infty$ for a certain $\mathit{\delta }$. Numerically, ${\mathit{\lambda }}_n$ would be even less than the computing precision when $n$ is large, which means that even without noise in measurements, the result would become inaccurate.

By our truncated SVD method, an appropriate $\mathrm{\Lambda }$ could be chosen so that $\frac{\mathit{\delta }}{\mathrm{\Lambda }}$ is bounded. In practice, we choose $\mathrm{\Lambda }={\mathit{\delta }}^{\mathit{\alpha }},$$0<\mathit{\alpha }<1$, then $\frac{\mathit{\delta }}{\mathrm{\Lambda }}+\mathrm{\Lambda }|z|\sim {\mathit{\delta }}^{1-\mathit{\alpha }}+\tilde{c}{\mathit{\delta }}^{\mathit{\alpha }}$ is independent of $n$.

Here, it is obvious that $\Vert {\tilde{\widehat{q}}}_n^{\mathit{\delta }}-\widehat{q}{\Vert }_{L^2(D)}=2\mathit{\pi }{\Vert {\tilde{q}}_n^{\mathit{\delta }}-q\Vert }_{L^2(D)}$ is nothing but $2\mathit{\pi }\Vert K_1^{+}{\mathit{\varphi }}^{\mathit{\delta }}{-q\Vert }_{L^2(D)}$ in (2.192.19 $$\begin{array}{c}\hfill {∥q,-,K_1^{+},\mathit{\varphi }∥}_{L^2(D)}⩽C\frac{[({\mathit{\mu }}_2+{\mathit{\nu }}_2)\Vert K_1^{+}\mathit{\varphi }{\Vert }_2]{}^2}{1-({\mathit{\mu }}_2+{\mathit{\nu }}_2){\Vert K_1^{+}\mathit{\varphi }\Vert }_2}+\tilde{C}{\Vert (I-K_1^{+}{K}_1)q\Vert }_{L^2(D)},\end{array}$$2.19 ) with ${\mathit{\varphi }}^{\mathit{\delta }}$ the perturbed data of measurements $\mathit{\varphi }$. Thus we have obtained a more detailed error estimate of the numerical method of the linearized inverse problem.

4 Numerical experiments

In this section, we give some numerical examples to illustrate the effectiveness of moment method with truncated SVD. In the first subsection, we give an example that is simple but really difficult to handle by ordinary methods, for its severe ill-posedness. And in the second subsection, we apply the moment method with truncated SVD to the inverse scattering problem we discussed above.

4.1 Numerical example involving inverse Laplace transform

In this subsection, we give some numerical examples of a much simpler inverse Laplace transform problem, where we only know a few discretization data. Our purpose is to show the readers how ill-posed the problem is when evanescent waves are involved.

Figure 18. $2c=\mathit{\lambda }/8$, $b=\mathit{\lambda }/2$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-6}$.

Figure 19. $2c=\mathit{\lambda }/8$, $b=\mathit{\lambda }$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-6}$.

Consider4.1 $$\begin{array}{c}\hfill g(t)={\int }_0^1{e}^{-ts}f(s)ds,t\in [0,1],\end{array}$$4.1 where we should reconstruct $f(t),t\in [0,1]$ from $n$ observations of $g(t_i),i=1,\cdots ,n.$ Moreover, the observations $g(t_i)$’s are perturbed by noise, i.e. $g(t_i):=g(t_i)+\mathit{\sigma }\mathrm{r}and,i=1,\cdots ,n$, where rand gives uniformly distributed numbers in [0,1].

In the numerical experiment, we set $g(t)={\displaystyle \frac{1-e^{-t}}{t}}$, so that $f(t)=1$ is the true solution. And we set that $n=8,\mathit{\sigma }={10}^{-4}.$

Firstly, we show the direct discretization method and Tikhonov regularization of the direct discretization method in Figure 2. In the direct discretization method, we use the trapezoidal rule; and in Tikhonov regularization, we choose the regularization parameter to be ${10}^{-2}$. The solution of the direct discretization method has nothing to do with the true solution, while the Tikhonov regularization method seems better, but not satisfying, either.

Next, we use moment method to reconstruct $f(t)$. Figure 3(a) uses just moment method, which is also a regularization method, but the result is much worse than what we expect. Figure 3(b) uses moment method with truncated SVD, and we get a satisfactory result. The truncated number is chosen as $\mathrm{\Lambda }={10}^{-4}$. We will apply the moment method with truncated SVD to the inverse scattering problem in the next subsection.

4.2 Numerical examples of inverse scattering problem

In this subsection, we give the numerical examples of the reconstruction of the scatterer in the linearized problem.

The computing model contains the domain of $[-\mathit{\lambda },\mathit{\lambda }]\times [-\mathit{\lambda },\mathit{\lambda }]$, where $\mathit{\lambda }$ is the wavelength of the incident wave. The index of refraction $n(\mathit{x})=1$ in the upper half plane and $n_0$ in the lower one. Two scatterers are presented on the surface of the substrate, in the domain $[-d,-c]\times [0,d-c]$ and $[c,d]\times [0,d-c]$, respectively, where $d-c=\mathit{\lambda }/16,q(\mathit{x})=1/2$, as is shown in Figure 4. The measurements of scattering field $u^s$ are collected on the segment $\mathrm{\Gamma }:=\{\mathit{x}:x_2=b,-L<x_1<L\}$. Here, we set $L=\mathit{\lambda }$. And the measurements data are perturbed by a noise level $\mathit{\sigma }$, i.e.$$u^s{{|}_{\mathrm{\Gamma }}:=u^s|}_{\mathrm{\Gamma }}(1+\mathit{\sigma }\mathrm{r}and),$$where rand gives uniformly distributed numbers in $[0,1]$. Thus, $|\mathit{\beta }-{\mathit{\beta }}^{\mathit{\delta }}|<c\mathit{\sigma }=\mathit{\delta }.$

Note that though we may set $\mathit{\sigma }=0$ in some experiments, it does not mean that there is no noise, but that the machine precision plays the role of noise up to $2\times {10}^{-16}$ which is not negligible in severely ill-posed problems.

Figure 20. $n_0=1.5,$$2c=\mathit{\lambda }/8$, $b=\mathit{\lambda }/2$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-6}$.

Figure 21. $n_0=1.5,$$2c=\mathit{\lambda }/8$, $b=\mathit{\lambda }$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-6}$.

Figure 22. $\mathit{\theta }\in [-\mathit{\pi }/2,-{\mathit{\theta }}_{cr}]\cup [{\mathit{\theta }}_{cr},\mathit{\pi }/2],$$2c=\mathit{\lambda }/8$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-6}$.

Example 1

In this example, we mainly show how to choose the truncated number $\mathrm{\Lambda }$ depending on the noise level $\mathit{\sigma }$, and investigate the influence of the measurements height $b$. In this example, we set $n_0=1.2$ and use scattering data by incident angles from $-\mathit{\pi }/2$ to $-\mathit{\pi }/2$, and do several experiments to reconstruct the scatterer by scattering field from $b=\mathit{\lambda }/2$ and $b=\mathit{\lambda }$, with the scatterers being apart by $2c=\frac{\mathit{\lambda }}{2},\frac{\mathit{\lambda }}{4},\frac{\mathit{\lambda }}{8}$, respectively.

Figures 5–10 show where the two scatterers lie apart by $\frac{\mathit{\lambda }}{2}$, and the scattering field is detected on the segment $\mathrm{\Gamma }$ with $b=\frac{\mathit{\lambda }}{2}$. Figures 5 and 6 show the results when $\mathit{\sigma }=0$, and taking the truncated number $\mathrm{\Lambda }={10}^{-8}$ and ${10}^{-6}$, respectively. It is readily seen that a truncated number should be small when noise level is low.

Next, we turn to the examples when the noise level increases. Figures 7 and 8 show the reconstruction result when the noise level is $\mathit{\sigma }={10}^{-6}$, and the truncated numbers are $\mathrm{\Lambda }={10}^{-6}$ and $\mathrm{\Lambda }={10}^{-3}$, respectively.

Figures 9 and 10 show the reconstruction results when the noise level is $\mathit{\sigma }={10}^{-3}$, and the truncated numbers are $\mathrm{\Lambda }={10}^{-6}$ and $\mathrm{\Lambda }={10}^{-3}$, respectively. It is clear in Figure 9 that the noise fails the reconstruction of the scatterer, but an appropriate truncated number in Figure 10 retrieves the resolution.

Now we continue the experiments when $b=\mathit{\lambda }$, with all the other data the same as before. Figures 11–13 show $(\mathit{\sigma },\mathrm{\Lambda })=(0,{10}^{-6})$, $(\mathit{\sigma },\mathrm{\Lambda })=({10}^{-6},{10}^{-6})$ and $(\mathit{\sigma },\mathrm{\Lambda })=({10}^{-3},{10}^{-3})$, respectively. Compared to those when $b=\mathit{\lambda }/2$, the reconstructions seem to be poor, but they can also give an acceptable resolution.

Next, we reduce the distance between the two scatterers to $\mathit{\lambda }/4$. Figures 14 and 15 show the results when noise levels are $\mathit{\sigma }=0$ and $\mathit{\sigma }={10}^{-3}$, with measurements height $b=\mathit{\lambda }/2$. Figures 16 and 17 show the results when noise levels are $\mathit{\sigma }=0$ and $\mathit{\sigma }={10}^{-3}$, with measurements height $b=\mathit{\lambda }$.

In the end of this example, when the distance between the two scatterers is reduced to $\mathit{\lambda }/8$, the reconstruction results only show the relative position of the scatterers but unable to clearly distinguish them, as is shown in Figures 18 and 19.

Example 2

In this example, we show how the index of fraction in the lower half-space $n_0$ influences the reconstruction. Here, the distance between the two scatters is still $\mathit{\lambda }/8$, which cannot distinguish the scatterers in the previous example. But when we set $n_0=1.5$ instead of 1.2, the reconstruction is much better. The results are shown in Figures 20 and 21, with measurements height $\mathit{\lambda }/2$ and $\mathit{\lambda }$, respectively. Especially, when the measurements height is $\mathit{\lambda }/2$, the two scatters can be distinguished clearly.

Figure 23. $\mathit{\theta }\in [-{\mathit{\theta }}_{cr},{\mathit{\theta }}_{cr}]$, $2c=\mathit{\lambda }/8$, $\mathit{\sigma }=0$, $\mathrm{\Lambda }={10}^{-6}$.

Example 3

In this example, we show how the incident angle influences the reconstruction. Here, we use all the parameters in Example 2 when the measurements height is $\mathit{\lambda }/2$. The critical angle when $u^t$ becomes evanescent wave is ${\mathit{\theta }}_{cr}=arcsin(1/n_0)$. So, we test the reconstruction when the incident angle varies between $[-\mathit{\pi }/2,-{\mathit{\theta }}_{cr}]\cup [{\mathit{\theta }}_{cr},\mathit{\pi }/2]$ in Figure 22, and when the incident angle varies between $[-{\mathit{\theta }}_{cr},{\mathit{\theta }}_{cr}]$ in Figure 23. The results tell us that fewer incident angles make the reconstruction poorer, but if we only use the evanescent waves, the results would also be acceptable.

The experiments show that the close measurements can lead to better resolution, and the truncated number plays a very important role in the reconstruction. In practice, we can get a good resolution when the truncated number is chosen close to the noise level. Besides, the larger refraction index of the lower half-space can improve the resolution of the reconstruction. More incident angles make the reconstruction better, but if we only use the evanescent waves, the results would also be acceptable for small scatterers. So, we can conclude that the evanescent waves play a very important role in near-field optics. But the exponentially ill-posedness generated by evanescent waves should be handled carefully by special regularization method, for example, the truncated SVD combined with moment method in this paper.

5 Conclusion

The inverse scattering problem in two-layered background is common and important in near-field optics. Based on the model of PSTM, we give an error bound of the inverse problem. We discuss inverse linear problem in detail and give the moment method with truncated SVD to handle the severely ill-posed problem with noisy data. Numerical experiments show that the method can give a resolution much better than $\mathit{\lambda }/2$.

There are two important future directions of the present study. The first concerns the more complicated three-dimensional model, where the three-dimensional Fredholm integral equation of the first kind should be solved. Another research direction is concerned with a finer error estimate of the Born approximation theoretically in our model.

Acknowledgments

The authors would like to thank the referees for their careful reading and valuable comments which led to improvement of the quality of our manuscript.

Notes

1 The research was supported by the National Natural Science Foundation of China (NSFC) [grant number 10971083], [grant number 11371172].