A stable sampling method for inverse scattering from periodic media

1 Introduction

In this paper we aim to numerically solve the inverse scattering problem for periodic media in $R^{2}$ . The periodic media of interest are unboundedly periodic in the horizontal direction and bounded in the vertical direction. These periodic media are motivated by one-dimensional photonic crystals and the inverse problem of interest is inspired by applications of nondestructive evaluations for photonic crystals. There have been an increasing amount of studies on numerical methods for shape reconstruction of periodic scattering media during the past years, see [Arens2005, Arens2003a, Bao2014, Elsch2012, Hadda2017, Lechl2013b, Nguye2014, Sandf2010, Yang2012, Nguye2016, Jiang2017, Cakon2019, Nguye2022b, Nguye2022]. Two major approaches that were studied in these papers are the factorization method and the near field imaging method. The latter method, which relies on a transformed field expansion, can provide super-resolved resolution. However, this method requires the periodic scattering structure to be a smooth periodic function multiplied by a small surface deformation parameter. The factorization method can essentially work for periodic scattering structures of arbitrary shape but it is not very robust against the noise in the scattering data. This method belongs to the class of sampling or qualitative methods that were introduced by D. Colton and A. Kirsch [Colto1996, Kirsc1998]. The factorization method aims to construct a necessary and sufficient characterization of the unknown scatterer from multi-static data. We refer to [Kirsc2008] for more details about the factorization method.

In this work we develop a sampling method with a new imaging functional to solve the inverse scattering problem for periodic media. This sampling method is inspired by the orthogonality sampling method [Potth2010] and the direct sampling method [Ito2012]. These two sampling methods share similar ideas and features and were studied independently. For simplicity we will refer to them as the orthogonality sampling method. The computation of the new imaging functional is simple and fast as one only needs to evaluate two finite sums involving the propagating modes of the scattered field data. Like the orthogonality sampling method, the proposed sampling method is also regularization free and very robust against noise in the data. The resolution of the new imaging functional is studied using Green’s identities and the Rayleigh expansion of the $α$ -quasiperiodic fields of the scattering problem. The stability of the functional is also established. The performance of the new imaging functional is studied in various contexts in the numerical study. The numerical study also shows that the proposed sampling method is more robust than the factorization method and more efficient than the orthogonality sampling method in reconstructing periodic scattering media. We also want to mention that although the orthogonality sampling method has been studied for inverse scattering from bounded objects [Potth2010, Gries2011, Ito2012, Ito2013, Kang2018, Harri2020, Le2022], its application to the inverse scattering problem for periodic media is still not known.

The paper is organized as follows. The basics of the scattering from periodic media and the inverse problem of interest are described in Section 2. The new imaging functional and its resolution and stability analysis are discussed in Section 3. Section 4 is dedicated to a numerical study of the new imaging functional and its comparison to the factorization method and the orthogonality sampling method.

2 Problem setup

We consider a two-dimensional medium which is unboundedly $2 π$ -periodic in $x_{1}$ -direction and bounded in $x_{2}$ -direction. Let $n$ be a bounded function which is $2 π$ -periodic with respect to $x_{1}$ . Suppose that the interior of the periodic medium is characterized by $n$ and that the exterior of the periodic medium is homogeneous which means $n = 1$ in these areas. Note that the period can be any arbitrary value, but it is chosen to be $2 π$ for the convenience of the presentation. For $α \in R$ , we define that a function $f$ is $α$ -quasiperiodic in $x_{1}$ if

From now on we will call functions with this property $α$ -quasiperiodic functions for short. A typical example of $α$ -quasiperiodic functions is a plane wave (e.g. $exp (i k (d_{1} x_{1} + d_{2} x_{2}))$ with $d_{1}^{2} + d_{2}^{2} = 1$ , $k > 0$ ). Suppose that the periodic medium is illuminated by an $α$ -quasiperiodic incident field $u_{i n}$ with wave number $k > 0$ . Note that since the medium is unboundedly periodic in $x_{1}$ , we are only interested incident fields propagating downward or upward toward the medium. The scattering of this incident field by the periodic medium produces the scattered field $u_{s c}$ described by

Δ u_{s c} + k^{2} u_{s c} = - k^{2} q u in% R^{2},

(1)

where $q$ is the contrast given by

q = n - 1.

It is well known for this scattering problem that the scattered field $u_{s c}$ must also be $α$ -quasiperiodic, and that the direct problem of finding the scattered field can be reduced to one period

Ω := (- π, π) \times R .

Let $D = supp (q) \cap Ω$ . For $h > 0$ such that

h > sup {| x_{2} | : (x_{1}, x_{2})^{⊤} \in D},

(2)

the direct scattering problem is completed by the Rayleigh expansion condition for the scattered field

u_{s c} (x) = ⎧ ⎨ ⎩ \begin{matrix} \sum_{j \in Z} u_{j}^{+} e^{i α_{j} x_{1} + i β_{j} (x_{2} - h)}, & x_{2} \geq h, \sum_{j \in Z} u_{j}^{-} e^{i α_{j} x_{1} - i β_{j} (x_{2} + h)}, & x_{2} \leq - h, \end{matrix}

(3)

where

α_{j} := α + j, β_{j} := ⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} \sqrt{k^{2} - α_{j}^{2}}, & k^{2} \geq α_{j}^{2} i \sqrt{α_{j}^{2} - k^{2}}, & k^{2} < α_{j}^{2} \end{matrix}, j \in Z,

and $(u_{j}^{\pm})_{j \in Z}$ are the (complex-valued) Rayleigh sequences of the scattered field $u_{s c}$ . The condition (3) means that the scattered field $u_{s c}$ is an outgoing wave. Note that only a finite number of terms in (3) are propagating plane waves which are called propagating modes, the rest are evanescent modes which correspond to exponentially decaying terms. From now, we call a function satisfying (3) a radiating function. In addition, we also assume that $β_{j}$ is nonzero for all $j$ which means the Wood anomalies are excluded in our analysis.

Well-posedness of scattering problem (1)–(3) is well-known, see for instance [Bonne1994]. For $r > 0$ define

Ω_{r} := (- π, π) \times (- r, r), Γ_{\pm r} := (- π, π) \times {\pm r} .

Recall the constant $h$ in (2). For the inverse problem of interest we measure the scattered field $u_{s c}$ on $Γ_{\pm r}$ for some $r \geq h$ . From the Rayleigh expansion of $u_{s c}$ , knowing $u_{s c}$ on $Γ_{\pm r}$ is equivalent to knowing the Rayleigh coefficients $(u_{j}^{\pm})_{j \in Z}$ . Since the evanescent modes are associated with the exponentially decaying terms in the Rayleigh expansion (3), it is typically difficult to obtain these modes in practice unless one can measure extremely near (e.g. within one wavelength) the periodic scatterers. Thus, we consider only the propagating modes for the scattering data of our inverse problem as follows.

Inverse problem. Given the Rayleigh coefficients $(u_{j}^{\pm})$ for $j \in Z$ such that $β_{j} > 0$ , determine $D$ .

3 The imaging functional and its properties

Recall that the $α$ -quasiperiodic Green function of the direct problem is given by

G (x, y) = \frac{i}{4 π} \sum j \in Z \frac{1}{β_{j}} e^{i α_{j} (x_{1} - y_{1}) + i β_{j} | x_{2} - y_{2} |}, x, y \in Ω, x_{2} \neq y_{2} .

(4)

It is well known that the direct problem is equivalent to the Lippmann-Schwinger equation

u_{s c} (x) = k^{2} \int_{D} G (x, y) q (y) u (y) d y .

(5)

Let $N$ be the number of incident fields we use for the inverse problem. Let $u_{j}^{\pm} (l)$ be the Rayleigh coefficients of the scattered field $u_{s c} (\cdot, l)$ generated by incident field $u_{i n} (\cdot, l)$ for $l = 1, 2, \dots, N$ . For $p \in N, z \in Ω$ , define the following imaging functional

where

g_{j}^{\pm} (z) = \frac{i}{4 π β_{j}} e^{- i α_{j} z_{1} \mp i β_{j} (z_{2} \mp h)} .

This functional $I (z)$ aims to determine $D$ in $Ω$ and $z$ plays the role of sampling points. We note that $g_{j}^{\pm} (z)$ are also the Rayleigh coefficients of the $α$ -quasiperiodic Green function. In the proof of following theorem we will drop the dependence of $u_{j}^{\pm}$ on $l$ for the convenience of the presentation but we keep $l$ in the total field $u (y, l)$ that is generated by incident field $u_{i n} (y, l)$ . The following theorem is the main theorem justifying the behavior of $I (z)$ as an imaging functional.

Theorem 1.

For $z \in Ω$ , the imaging functional satisfies

I (z) = {(\frac{k^{2}}{8 π})}^{p} N \sum l = 1 {∣ ∣ ∣ \int_{D} J_{0} (k | y - z |) q (y) u (y, l) d y + \int_{D} w (y, z) q (y) u (y, l) d y ∣ ∣ ∣}_{D}^{p}

where

	$w (y, z) = \sum j \in Z ∖ {0} cos (α j 2 π) J_{0} (k \sqrt{(y_{1} - z_{1} + 2 j π)^{2} + (y_{2} - z_{2})^{2}})$
	$- \sum j \in Z ∖ {0} sin (α j 2 π) Y_{0} (k \sqrt{(y_{1} - z_{1} + 2 j π)^{2} + (y_{2} - z_{2})^{2}})$		(6)

where $J_{0}$ and $Y_{0}$ are respectively the Bessel functions of first and second kinds.

Remark 2.

The series in (1) converge since $y_{1} - z_{1} + 2 j π$ is always nonzero for $y_{1}, z_{1} \in (- π, π)$ and $j \neq 0$ . We know that $J_{0} (t)$ and $Y_{0} (t)$ decay like $O (1 / \sqrt{t})$ as $t \to \infty$ . Thus the terms with larger $| j |$ in (1) are small terms. The terms with smaller $| j |$ can be also small with a suitably larger wave number $k$ (e.g. the choice of $k = 2 π$ works well in our numerical study). Therefore, the integral associated with $w (y, z)$ is expected to be small compared with the integral with the kernel $J_{0}$ . Thanks to the behavior of $J_{0}$ the latter integral is expected to have small values as $z$ is outside $D$ and to have much larger values as $z$ is inside $D$ .

Proof.

For any $x_{s}, x_{t} \in Ω_{h}$ , $G (x, x_{s})$ solves

Δ G (x, x_{s}) + k^{2} G (x, x_{s}) = - δ (x - x_{s}), x \in Ω_{h} .

Multiplying both sides by $¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t})$ and integrating over $Ω_{h}$ gives

\int_{Ω_{h}} (Δ G (x, x_{s}) + k^{2} G (x, x_{s})) ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) d x = - ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x_{s}, x_{t}) .

(7)

Similarly, $¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t})$ solves

Δ ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) + k^{2} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) = - δ (x - x_{t}), x \in Ω_{h},

thus by multiplying both sides by $G (x, x_{s})$ and integrating over $Ω_{h}$ we obtain

\int_{Ω_{h}} (Δ ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) + k^{2} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t})) G (x, x_{s}) d x = - G (x_{s}, x_{t}) .

(8)

Subtracting (8) from (7) yields

	$G (x_{s}, x_{t}) - ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x_{s}, x_{t})$	$= \int_{Ω_{h}} Δ G (x, x_{s}) ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) - Δ ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) G (x, x_{s}) d x$
		$= \int_{\partial Ω_{h}} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) \frac{\partial G (x, x_{s})}{\partial n} - G (x, x_{s}) \frac{\partial ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t})}{\partial n} d s (x),$

where $n$

is an unit normal outward vector to

\partial Ω_{h}

. Note that

\partial Ω_{h} = Γ \cup {\pm π} \times (- h, h)

, but thanks to the

α

-quasiperiodicity of

G

\int_{{\pm π} \times (- h, h)} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) \frac{\partial G (x, x_{s})}{\partial n} - G (x, x_{s}) \frac{\partial ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t})}{\partial n} d s (x) = 0.

Therefore

G (x_{s}, x_{t}) - ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x_{s}, x_{t}) = \int_{Γ_{h} \cup Γ_{- h}} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) \frac{\partial G (x, x_{s})}{\partial n} - G (x, x_{s}) \frac{\partial ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t})}{\partial n} d s (x) .

(9)

The left-hand side of (9) equals $2 i I m G (x_{s}, x_{t})$ , now we calculate the right-hand side. Since $G$ satisfies the radiation condition, for a fixed $y \in Ω_{h}$ ,

G (x, y) = \sum j \in Z g_{j}^{\pm} (y) e^{i α_{j} x_{1}}, x \in Γ_{\pm h}

(10)

and

\frac{\partial G (x, x_{t})}{\partial n} = \sum j \in Z i β_{j} g_{j}^{\pm} (x_{s}) e^{i α_{j} x_{1}}, x \in Γ_{\pm h} .

Therefore

	$\int_{Γ_{\pm h}} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) \frac{\partial G (x, x_{s})}{\partial n} d s (x)$	$= \int_{- π}^{π} \sum j_{1}, j_{2} \in Z i β_{j_{2}} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j_{1}}^{\pm} (x_{t}) g_{j_{2}}^{\pm} (x_{s}) e^{i (α_{j_{2}} - α_{j_{1}}) x_{1}} d x_{1}$
		$= \sum j_{1}, j_{2} \in Z i β_{j_{2}} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j_{1}}^{\pm} (x_{t}) g_{j_{2}}^{\pm} (x_{s}) \int_{- π}^{π} e^{i (j_{2} - j_{1}) x_{1}} d x_{1},$

and since

\int_{- π}^{π} e^{i (j_{2} - j_{1}) x_{1}} d x_{1} = {\begin{matrix} 2 π, & j_{1} = j_{2} 0, & j_{1} \neq j_{2} \end{matrix}

we have

\int_{Γ_{h} \cup Γ_{- h}} ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t}) \frac{\partial G (x, x_{s})}{\partial n} d s (x) = 2 π i \sum j \in Z β_{j} (¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (x_{t}) g_{j}^{+} (x_{s}) + ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (x_{t}) g_{j}^{-} (x_{s})) .

(11)

Similarly

\int_{Γ_{h} \cup Γ_{- h}} G (x, x_{s}) \frac{\partial ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, x_{t})}{} \partial n d s (x) = - 2 π i \sum j \in Z_{j} (¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (x_{t}) g_{j}^{+} (x_{s}) + ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (x_{t}) g_{j}^{-} (x_{s})) .

(12)

Combining (11) and (12) yields

G (x_{s}, x_{t}) - ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x_{s}, x_{t}) = 2 π i \sum j \in Z (β_{j} +_{j}) (¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (x_{t}) g_{j}^{+} (x_{s}) + ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (x_{t}) g_{j}^{-} (x_{s})),

that is

I m G (x_{s}, x_{t}) = 2 π \sum j \in Z R e β_{j} (¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (x_{t}) g_{j}^{+} (x_{s}) + ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (x_{t}) g_{j}^{-} (x_{s})) .

(13)

From the Rayleigh expansion (3) we have that

u_{j}^{\pm} = \frac{1}{2 π} \int_{Γ_{\pm h}} u_{s c} (x) e^{- i α_{j} x_{1}} d s (x) .

Thus, for $j \in Z$ , substituting (5) into the integral above gives

	$u_{j}^{\pm}$	$= \frac{1}{2 π} \int_{Γ_{\pm h}} (k^{2} \int_{D} G (x, y) q (y) u (y) d y) e^{- i α_{j} x_{1}} d s (x)$
		$= k^{2} \int_{D} (\frac{1}{2 π} \int_{Γ_{\pm}} G (x, y) e^{- i α_{j} x_{1}} d s (x)) q (y) u (y) d y = k^{2} \int_{D} g_{j}^{\pm} (y) q (y) u (y) d y .$

For $z \in Ω$ , using the Rayleigh coefficients of $u_{j}^{\pm}$ above and (13) we obtain

	$2 π \sum j \in Z R e β_{j} (u_{j}^{+} ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (z) + u_{j}^{-} ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (z))$
	$= k^{2} \int_{D} 2 π \sum j \in Z R e β_{j} (g_{j}^{+} (y) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (z) + g_{j}^{-} (y) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (z)) q (y) u (y) d y$
	$= k^{2} \int_{D} I m G (y, z) q (y) u (y) d y .$

Now the proof follows from using another representation of the $α$ -quasiperiodic Green function

G (y, z) = \sum j \in Z e^{- i α_{j} 2 π} \frac{i}{4} H_{0}^{(1)} (k \sqrt{(y_{1} - z_{1} + 2 j π)^{2} + (y_{2} - z_{2})^{2}})

where the series converges for $(y_{1} - z_{1}, y_{2} - z_{2}) \neq (2 j π, 0)$ , $j \in Z$ . ∎

In the next theorem we will establish a stability estimate for the imaging functional. Assume

u_{s c} (\cdot, l) \in L^{2} (Γ_{r} \cup Γ_{- r})

for all

l = 1, \dots, N

Theorem 3.

For $δ > 0$ , denote by $u_{s c, δ}$ and ${(u_{δ, j}^{\pm})}_{j \in Z}$ the noisy scattered wave and its Rayleigh sequences respectively, for which we have

N \sum l = 1 ∥ u_{s c, δ} (\cdot, l) - u_{s c} (\cdot, l) ∥_{L^{2} (Γ_{r} \cup Γ_{- r})} \leq δ .

Define

I_{δ} (z) := N \sum l = 1 {∣ ∣ ∣ ∣ \sum j : β_{j} > 0 β_{j} (u_{δ, j}^{+} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (z) + u_{δ, j}^{-} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (z)) ∣ ∣ ∣ ∣}_{j}^{p} .

Then the following stability property holds

| I_{δ} (z) - I (z) | = O (δ), as δ \to 0

for every $z \in Ω$ .

Proof.

For $l = 1, \dots, N$ and $j \in Z$ such that $β_{j} > 0$ , we have

| u_{δ, j}^{\pm} (l) - u_{j}^{\pm} (l) | \leq \frac{1}{2 π} \int_{Γ_{\pm r}} | u_{s c}^{δ} (x, l) - u_{s c} (x, l) | d s (x),

thus, by Cauchy-Schwarz inequality,

N \sum l = 1 | u_{δ, j}^{\pm} (l) - u_{j}^{\pm} (l) | \leq N \sum l = 1 ∥ u_{s c}^{δ} (\cdot, l) - u_{s c} (\cdot, l) ∥_{L^{2} (Γ_{r} \cup Γ_{- r})} \leq δ .

Therefore

	$∣ ∣ ∣ ∣ N \sum l = 1 ∣ ∣ ∣ ∣ \sum j : β_{j} > 0 β_{j} (u_{δ, j}^{+} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (z) + u_{δ, j}^{-} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (z)) ∣ ∣ ∣ ∣ - N \sum l = 1 ∣ ∣ ∣ ∣ \sum j : β_{j} > 0 β_{j} (u_{j}^{+} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (z) + u_{j}^{-} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (z)) ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣$
	$\leq N \sum l = 1 ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ \sum j : β_{j} > 0 β_{j} (u_{δ, j}^{+} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (z) + u_{δ, j}^{-} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (z)) ∣ ∣ ∣ ∣ - ∣ ∣ ∣ ∣ \sum j : β_{j} > 0 β_{j} (u_{j}^{+} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (z) + u_{j}^{-} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (z)) ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣$
	$\leq N \sum l = 1 \sum j : β_{j} > 0 β_{j} (∣ ∣ u_{δ, j}^{+} (l) - u_{j}^{+} (l) ∣ ∣ ∣ ∣ g_{j}^{+} (z) ∣ ∣ + ∣ ∣ u_{δ, j}^{-} (l) - u_{j}^{-} (l) ∣ ∣ ∣ ∣ g_{j}^{-} (z) ∣ ∣)$
	$\leq C δ,$

where

C = (max j : β_{j} > 0 {∥ ∥ g_{j}^{+} ∥ ∥}_{L^{\infty} (Ω)} + max j : β_{j} > 0 {∥ ∥ g_{j}^{-} ∥ ∥}_{L^{\infty} (Ω)}) \sum j : β_{j} > 0 β_{j} .

For brevity, set

	$a_{l}$
	$b_{l}$	$:= ∣ ∣ ∣ ∣ \sum j : β_{j} > 0 β_{j} (u_{j}^{+} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (z) + u_{j}^{-} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (z)) ∣ ∣ ∣ ∣,$
	$γ$	$:= max \begin{matrix} j : β_{j} > 0 l = 1, \dots, N \end{matrix} {∣ ∣ u_{j}^{+} (l) ∣ ∣, ∣ ∣ u_{j}^{-} (l) ∣ ∣},$

then we have

N \sum l = 1 | a_{l} - b_{l} | \leq C δ, b_{l} \leq C γ,

for all $l = 1, \dots, N$ . Hence, for all $z \in Ω$ ,

	$\| I_{δ} (z) - I (z) \|$	$\leq N \sum l = 1 ∣ ∣ a_{l}^{p} - b_{l}^{p} ∣ ∣ = N \sum l = 1 (\| a_{l} - b_{l} \| ∣ ∣ ∣ ∣ p - 1 \sum m = 0 a_{l}^{m} b_{l}^{p - 1 - m} ∣ ∣ ∣ ∣)$
		$\leq (N \sum l = 1 \| a_{l} - b_{l} \|) (N \sum l = 1 ∣ ∣ ∣ ∣ p - 1 \sum m = 0 a_{l}^{m} b_{l}^{p - 1 - m} ∣ ∣ ∣ ∣)$
		$\leq C δ N \sum l = 1 p - 1 \sum m = 0 (\| a_{l} - b_{l} \| + b_{l})^{m} b_{l}^{p - 1 - m}$
		$\leq C δ N \sum l = 1 p - 1 \sum m = 0 2^{m} ({\| a_{l} - b_{l} \|}_{l}^{m} + b_{l}^{m}) b_{l}^{p - 1 - m}$
		$\leq C δ p - 1 \sum m = 0 N \sum l = 1 2^{m} ({\| a_{l} - b_{l} \|}_{l}^{m} + γ^{m}) γ^{p - 1 - m}$
		$\leq C δ p - 1 \sum m = 0 2^{m} (C^{m} δ^{m} γ^{p - 1 - m} + N γ^{p - 1})$
		$= O (δ), as δ \to 0.$

This completes the proof. ∎

4 Numerical study

In this section, we test the performance of the proposed sampling method with respect to the number of incident sources, the levels of noise in the data, the wave numbers, the values of parameter $α$ (in (3)) and exponent $p$ , and the shape of the periodic scatterers. For the latter category we compare the performance of the proposed sampling method with those of the factorization method and the orthogonality sampling method.

In all of the numerical examples below we choose the following parameters:

	$sampling domain = (- π, π) \times (- 1, 1),$
	$measurement boundary Γ_{\pm 2} = (- π, π) \times {\pm 2},$
	$location of incident sources Γ_{\pm 3} = (- π, π) \times {\pm 3} .$

The sampling domain is uniformly discretized in each dimension with $128 \times 96$ sampling points. The boundary measurements $Γ_{\pm 2}$ are discretized uniformly with 64 points on each boundary. By using $N$ incident sources we mean to consider

u_{i n} (x, l) = G (x, s_{l}), x \in Ω, s_{l} \in Γ_{\pm 3}, l = 1, \dots, N,

where $N / 2$ sources are uniformly located on $Γ_{+ 3}$ and $N / 2$ sources are uniformly located on $Γ_{- 3}$ . We generate the synthetic scattering data by solving the direct problem with the spectral Galerkin method studied in [Lechl2014]. With artificial noise added to the synthetic scattering data we implement the imaging functional

I_{δ} (z) = N \sum l = 1 {∣ ∣ ∣ ∣ \sum j : β_{j} > 0 β_{j} (u_{δ, j}^{+} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{+} (z) + u_{δ, j}^{-} (l) ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ g_{j}^{-} (z)) ∣ ∣ ∣ ∣}_{j}^{p} .

To compare with the orthogonality sampling method we implement following imaging functional

I_{O S M} (z) = N \sum l = 1 {∣ ∣ ∣ \int_{Γ_{+ 2} \cup Γ_{- 2}} u_{s c, δ} (x) ¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ G (x, z) d s (x) ∣ ∣ ∣}_{Γ_{+ 2} \cup Γ_{- 2}}^{p} .

This imaging functional of the orthogonality sampling method can be rewritten in the modal form by using the Rayleigh expansion of the scattered field. In this modal form the evanescent modes or the exponentially decaying terms can be neglected and there remains only a finite number of the propagating modes. If we drop $β_{j}$ in $I_{δ} (z)$ , we will approximately obtain the modal form of $I_{O S M} (z)$ . We refer to [Nguye2014] for the imaging functional of the factorization method implemented for the test in this section.

4.1 Reconstruction with one incident source (Figure 1)

In this section we test the performance of the sampling method for data generated by only one incident source. Here the parameters are chosen as $α = 0$ , wave number $k = 2 π$ , and exponent $p = 4$ . We add 20 $%$ artificial noise to the scattered field data ( $δ = 20 %$ ). From Figure 1 we can see that the method is able to reconstruct small scatterers quite well. This capability is an advantage over classical sampling methods (e.g. linear sampling method, factorization method) in terms of computational efficiency since classical methods are only able to reconstruct targets with data generated by multiple incident fields. However, the proposed sampling method fails to reconstruct scatterers with extended shape, which is reasonable since only one incident source and one wave number are used for the data. We refer to Figure 2 for improved results when multiple incident sources are used for the reconstruction.

4.2 Reconstruction with multiple incident sources (Figure 2)

In this section we test the performance of the method for reconstructing extended scatterers with different number of incident sources. The parameters are the same as in Figure 1, meaning $α = 0$ , wave number $k = 2 π$ , exponent $p = 4$ , and 20 $%$ artificial noise is added to the scattered field data. If in Figure 1 the method fails to reconstruct extended ellipses with one incident source, the reconstructions are improved with more sources, see Figure 2. The reconstructions already look reasonable with 32 incident sources and continue to improve when more sources are used. The results remain almost the same even if more than 128 incident sources are used to generate the scattering data.

4.3 Reconstruction with different levels of noise in the data (Figure 3)

In this section we test the performance of the method for reconstructing extended scatterers with different levels of artificial noise in the data ( $10 %, 20 %$ and $40 %$ ). Here the parameters are chosen as $α = 0$ , wave number $k = 2 π$ , exponent $p = 4$ , and 128 incident sources are used to generate the scattering data. It can be seen from Figure 3 that all reconstructions are not affected by the amounts of noise added to the data. Furthermore, the reconstructions will remain essentially the same even with much higher amounts of noise in the data. This is not a surprise since the evanescent modes are typically sensitive with noise but the sampling method only uses propagating modes. The great robustness against noise in the data was also seen in the orthogonality sampling methods studied in [Potth2010, Harri2020, Le2022].

4.4 Reconstruction with different wave numbers (Figure 4)

In this section we test the performance of the method for reconstructing extended scatterers with different wave numbers ( $k = π, 2 π, 3 π$ ). Here the parameters are chosen as $α = 0$ , exponent $p = 4$ , 128 incident sources are used to generate the scattering data, and $20 %$ noise added to the data. We can see from Figure 4 that the resolution of reconstructions improve as the wave number increases, which is reasonable. However, it is also known that as the wave number increases, the inverse problem will become more difficult to deal with. We can notice some small effect of larger wave numbers in the reconstruction for $k = 3 π$ .

4.5 Reconstruction with different values of $α$ and $p$ (Figures 5–6)

In this section we test the performance of the method for different values of the parameters $α$ and $p$ . The other parameters are chosen as $k = 2 π$ , 128 incident sources are used to generate the scattering data, and $20 %$ noise added to the data. From Figures 5–6, we can see that the sampling method works well for different values of $α$ . Although the reconstructions look reasonable for exponents $p = 2$ and $p = 3$ , the exponent $p = 4$ seems to be an ideal exponent for the imaging functional. With larger $p$ the reconstructions will have cleaner background but lose some small details of the scatterers.

4.6 Reconstruction of different shapes and comparison with the factorization method and the orthogonality sampling method (Figures 7–10)

In this section we test the performance of the proposed sampling method for different shapes of periodic scatterers and compare its performance with those of the factorization method and the orthogonality sampling method. The parameters are chosen as $α = 0$ , wave number $k = 2 π$ , exponent $p = 4$ , 128 incident sources are used to generate the scattering data, and $20 %$ noise added to the data. It is obvious from all reconstructions in Figures 7–10 that the factorization method suffers severely from the $20 %$ amount of noise added to the scattering data. Actually, the reconstructions of the factorization method can be greatly affected even with smaller amounts of noises (e.g. $5 %$ or $7 %$ ). This unstable behavior of the factorization method in reconstructing periodic media was also reported in [Arens2005, Nguye2016]. The reconstructions of the orthogonality sampling method are as stable as those of the proposed sampling method but it is also clear from the pictures that the accuracy in the reconstructions of the orthogonality sampling method is much worse than that of the proposed sampling method. The proposed sampling method may provide reasonable reconstructions for different shapes considered in the test. This indicates a high efficiency of this sampling method.

Acknowledgement. This work was partially supported by NSF grant DMS-1812693.

	$\| I_{δ} (z) - I (z) \|$	$\leq N \sum l = 1 ∣ ∣ a_{l}^{p} - b_{l}^{p} ∣ ∣ = N \sum l = 1 (\| a_{l} - b_{l} \| ∣ ∣ ∣ ∣ p - 1 \sum m = 0 a_{l}^{m} b_{l}^{p - 1 - m} ∣ ∣ ∣ ∣)$
		$\leq (N \sum l = 1 \| a_{l} - b_{l} \|) (N \sum l = 1 ∣ ∣ ∣ ∣ p - 1 \sum m = 0 a_{l}^{m} b_{l}^{p - 1 - m} ∣ ∣ ∣ ∣)$
		$\leq C δ N \sum l = 1 p - 1 \sum m = 0 (\| a_{l} - b_{l} \| + b_{l})^{m} b_{l}^{p - 1 - m}$
		$\leq C δ N \sum l = 1 p - 1 \sum m = 0 2^{m} ({\| a_{l} - b_{l} \|}_{l}^{m} + b_{l}^{m}) b_{l}^{p - 1 - m}$
		$\leq C δ p - 1 \sum m = 0 N \sum l = 1 2^{m} ({\| a_{l} - b_{l} \|}_{l}^{m} + γ^{m}) γ^{p - 1 - m}$
		$\leq C δ p - 1 \sum m = 0 2^{m} (C^{m} δ^{m} γ^{p - 1 - m} + N γ^{p - 1})$
		$= O (δ), as δ \to 0.$

A stable sampling method for inverse scattering from periodic media

Authors

Imaging of bi-anisotropic periodic structures from electromagnetic near field data

Orthogonality sampling type methods for an inverse acoustic scattering problem

Inverse scattering by a random periodic structure

An inverse scattering approach for geometric body generation: a machine learning perspective

Monostatic sampling methods in limited-aperture configuration

Transfer Learning in General Lensless Imaging through Scattering Media

Data completion algorithms and their applications in inverse acoustic scattering with limited-aperture backscattering data

1 Introduction

2 Problem setup

3 The imaging functional and its properties

Theorem 1.

Remark 2.

Proof.

Theorem 3.

Proof.

4 Numerical study

4.1 Reconstruction with one incident source (Figure 1)

4.2 Reconstruction with multiple incident sources (Figure 2)

4.3 Reconstruction with different levels of noise in the data (Figure 3)

4.4 Reconstruction with different wave numbers (Figure 4)

4.5 Reconstruction with different values of $α$ and $p$ (Figures 5–6)

4.6 Reconstruction of different shapes and comparison with the factorization method and the orthogonality sampling method (Figures 7–10)

References

A stable sampling method for inverse scattering from periodic media

Authors

Related Research

Imaging of bi-anisotropic periodic structures from electromagnetic near field data

Orthogonality sampling type methods for an inverse acoustic scattering problem

Inverse scattering by a random periodic structure

An inverse scattering approach for geometric body generation: a machine learning perspective

Monostatic sampling methods in limited-aperture configuration

Transfer Learning in General Lensless Imaging through Scattering Media

Data completion algorithms and their applications in inverse acoustic scattering with limited-aperture backscattering data

This week in AI

1 Introduction

2 Problem setup

3 The imaging functional and its properties

Theorem 1.

Remark 2.

Proof.

Theorem 3.

Proof.

4 Numerical study

4.1 Reconstruction with one incident source (Figure 1)

4.2 Reconstruction with multiple incident sources (Figure 2)

4.3 Reconstruction with different levels of noise in the data (Figure 3)

4.4 Reconstruction with different wave numbers (Figure 4)

4.5 Reconstruction with different values of α and p (Figures 5–6)

4.6 Reconstruction of different shapes and comparison with the factorization method and the orthogonality sampling method (Figures 7–10)

References

4.5 Reconstruction with different values of $α$ and $p$ (Figures 5–6)