Near-field imaging of locally perturbed periodic surfaces

TL;DR

This paper develops a Floquet-Bloch transform-based numerical method to reconstruct both the periodic surface and local perturbations from near-field data, addressing challenges posed by non-quasi-periodic scattered fields.

Contribution

It introduces a novel numerical approach combining sampling and Newton-CG methods for inverse scattering of locally perturbed periodic surfaces.

Findings

Effective reconstruction of surface and perturbation demonstrated through numerical examples.

Method handles non-quasi-periodic scattered fields where classical methods fail.

Reconstruction accuracy depends on incident field properties.

Abstract

This paper concerns the inverse scattering problem to reconstruct a locally perturbed periodic surface. Different from scattering problems with quasi-periodic incident fields and periodic surfaces, the scattered fields are no longer quasi-periodic. Thus the classical method for quasi-periodic scattering problems no longer works. In this paper, we apply a Floquet-Bloch transform based numerical method to reconstruct both the unknown periodic part and the unknown local perturbation from the near-field data. By transforming the original scattering problem into one defined in an infinite rectangle, the information of the surface is included in the coefficients. The numerical scheme contains two steps. The first step is to obtain an initial guess, i.e., the locations of both the periodic surfaces and the local perturbations, from a sampling method. The second step is to reconstruct the surface. As is proved in this paper, for some incident fields, the corresponding scattered fields carry little information of the perturbation. In this case, we use this scattered field to reconstruct the periodic surface. Then we could apply the data that carries more information of the perturbation to reconstruct the local perturbation. The Newton-CG method is applied to solve the associated optimization problems. Numerical examples are given at the end of this paper to show the efficiency of the numerical method.

Equations

$$\Gamma :=\left\{(x_1,\zeta (x_1)):x_1\in \mathbb{R}\right\}.$$

$${\Gamma }_p:=\left\{(x_1,{\zeta }_p(x_1)):x_1\in \mathbb{R}\right\}.$$

$$\Delta u+k^{2}u=0\text{ in }{\Omega }_p.$$

$$u=0\text{ on }{\Gamma }_p.$$

$$u^s(x_1,x_2)=\frac{1}{2\pi }{\int }_{\mathbb{R}}{e}^{\mathrm{i}\xi x_1+\mathrm{i}\sqrt{k^2-{\xi }^2}(x_2-H)}{\hat{u}}^s(\xi ,H)\mathrm{d}\xi ,x_2\ge H,$$

$$\left(T^{+}\phi \right)(x_1)=\frac{\mathrm{i}}{2\pi }{\int }_{\mathbb{R}}\sqrt{k^2-{\xi }^2}{e}^{\mathrm{i}\xi x_1}\hat{\phi }(\xi )\mathrm{d}\xi ,\phi =\frac{1}{2\pi }{\int }_{\mathbb{R}}{e}^{\mathrm{i}\xi x_1}\hat{\phi }(\xi )\mathrm{d}\xi .$$

$$\frac{\partial u}{\partial x_{2}}(x_{1},H)=T^{+}\left[u\big{|}_{\Gamma_{H}}\right]+f\,\text{ on }\Gamma_{H},\,\text{where}\,\,f=\frac{\partial u^{i}}{\partial x_{2}}(x_{1},H)-T^{+}\left(u^{i}\big{|}_{\Gamma_{H}}\right).$$

$$H_r^1({\Omega }_H^p):=\left\{\phi \in {\mathcal{D}}^{\prime }({\Omega }_H^p):(1+|x{|}^2{)}^{r/2}\phi \in H^1({\Omega }_H^p)\right\}.$$

$$\int_{\Omega^{p}_{H}}\left[\nabla u\cdot\nabla\overline{v}-k^{2}u\overline{v}\right]\,\mathrm{d}x\,-\int_{\Gamma_{H}}T^{+}\left(u\big{|}_{\Gamma_{H}}\right)\overline{v}\,\mathrm{d}s\,=\int_{\Gamma_{H}}f\overline{v}\,\mathrm{d}s\,$$

$$W=(-\pi ,\pi ],W^{\ast }=(-1/2,1/2].$$

$${\mathcal{J}}_D\phi (\alpha ,x)=\sum _{j\in \mathbb{Z}}\phi \left(x+\left(\begin{array}{c}2\pi j\\ 0\end{array}\right)\right)e^{2\mathrm{i}\pi j\alpha }.$$

$$\Vert \phi {\Vert }_{H_0^r(W_{{\Lambda }^{\ast }};H_{\alpha }^s(D^{2\pi }))}={\left[\sum _{\gamma =0}^r{\int }_{W^{\ast }}\Vert {\partial }_{\alpha }^{\gamma }\psi (\alpha ,\cdot ){\Vert }_{H_{\alpha }^s(D^{2\pi })}^2\right]}^{1/2}.$$

$$({\mathcal{J}}_D^{-1}\phi )\left(x+\left(\begin{array}{c}2\pi j\\ 0\end{array}\right)\right)={\int }_{W^{\ast }}\phi (\alpha ,x)e^{2\mathrm{i}\pi j\alpha }\mathrm{d}\alpha ,x\in D^{2\pi },$$

$${\Phi }_{\zeta }:x\mapsto \left(x_1,x_2+\frac{(x_2-H_0{)}^3}{(h-H_0{)}^3}(\zeta (x_1)-h)\right);{\Phi }_{{\zeta }_p}:x\mapsto \left(x_1,x_2+\frac{(x_2-H_0{)}^3}{(h-H_0{)}^3}({\zeta }_p(x_1)-h)\right).$$

$$\int_{D}\left[A_{\zeta_{p}}\nabla u_{D}\cdot\nabla\overline{v_{D}}-k^{2}c_{\zeta_{p}}u_{D}\overline{v_{D}}\right]\,\mathrm{d}x\,-\int_{\Gamma_{H}}T^{+}\left(u_{D}\big{|}_{\Gamma_{H}}\right)\overline{v_{D}}\,\mathrm{d}s\,=\int_{\Gamma_{H}}f\overline{v_{D}}\,\mathrm{d}s\,,$$

$$A_{{\zeta }_p}(x)=\left|\det \nabla {\Phi }_{{\zeta }_p}(x)\right|\left[{\left(\nabla {\Phi }_{{\zeta }_p}(x)\right)}^{-1}{\left({\left(\nabla {\Phi }_{{\zeta }_p}(x)\right)}^{-1}\right)}^{\top }\right]\in L^{\infty }(D,{\mathbb{R}}^{2\times 2});$$

$$c_{{\zeta }_p}(x)=\left|\det \nabla {\Phi }_{{\zeta }_p}(x)\right|\in L^{\infty }(D).$$

$$A_{\zeta }(x)=\left|\det \nabla {\Phi }_{\zeta }(x)\right|\left[{\left(\nabla {\Phi }_{\zeta }(x)\right)}^{-1}{\left({\left(\nabla {\Phi }_{\zeta }(x)\right)}^{-1}\right)}^{\top }\right]\in L^{\infty }(D,{\mathbb{R}}^{2\times 2});$$

$$c_{\zeta }(x)=\left|\det \nabla {\Phi }_{\zeta }(x)\right|\in L^{\infty }(D).$$

$${\int }_{W^{\ast }}{a}_{\alpha }(w(\alpha ,\cdot ),z(\alpha ,\cdot ))\mathrm{d}\alpha +b({\mathcal{J}}_{\Omega }^{-1}w,{\mathcal{J}}_{\Omega }^{-1}z)={\int }_{W^{\ast }}{\int }_{{\Gamma }_H^{2\pi }}F(\alpha ,\cdot )\overline{z(\alpha ,\cdot )}\mathrm{d}\alpha \mathrm{d}s,$$

$$a_{\alpha }(u,v)={\int }_{{\Omega }_H^{2\pi }}\left[A_{\zeta }\nabla u\cdot \nabla \overline{v}-k^2{c}_{\zeta }u\overline{v}\right]\mathrm{d}x-{\int }_{{\Gamma }_H^{2\pi }}\left(T_{\alpha }^{+}u\right)\overline{v}\mathrm{d}s;$$

$$b(u,v)={\int }_{{\Omega }_H^{2\pi }}\left[(A_{{\zeta }_p}-A_{\zeta })\nabla u\cdot \nabla \overline{v}-k^2(c_{{\zeta }_p}-c_{\zeta })u\overline{v}\right]\mathrm{d}x;$$

$$F(\alpha ,\cdot )=\frac{\partial ({\mathcal{J}}_{\Omega }{u}^i)(\alpha ,\cdot )}{\partial x_2}-T_{\alpha }^{+}({\mathcal{J}}_{\Omega }{u}^i)(\alpha ,\cdot );$$

$$T_{\alpha }^{+}\phi =\mathrm{i}\sum _{j\in \mathbb{Z}}\sqrt{k^2-|j-\alpha {|}^2}\hat{\phi }(j)e^{\mathrm{i}(j-\alpha )x_1},\phi =\sum _{j\in \mathbb{Z}}\hat{\phi }(j)e^{\mathrm{i}(j-\alpha )x_1}.$$

$$\int_{\Omega_{H}}\left[\nabla u_{0}\cdot\nabla\overline{v}-k^{2}u_{0}\overline{v}\right]\,\mathrm{d}x\,-\int_{\Gamma_{H}}T^{+}\left(u_{0}\big{|}_{\Gamma_{H}}\right)\overline{v}\,\mathrm{d}s\,=\int_{\Gamma_{H}}f\overline{v}\,\mathrm{d}s\,.$$

$$\left|u_0(x_1,x_2)\right|\le C(1+x_1^2{)}^{-r/2-1/4}.$$

$$\int_{\Omega_{H}}\left[\nabla u_{0}^{L}\cdot\nabla\overline{v}-k^{2}u_{0}^{L}\overline{v}\right]\,\mathrm{d}x\,-\int_{\Gamma_{H}}T^{+}\left(u_{0}^{L}\big{|}_{\Gamma_{H}}\right)\overline{v}\,\mathrm{d}s\,=\int_{\Gamma_{H}}f_{L}\overline{v}\,\mathrm{d}s\,$$

$$\left|u_0(x_1+2\pi L,x_2)\right|\le C|2\pi L{|}^{-r-1/2},(x_1,x_2)\in {\Omega }_H^{2\pi },$$

$$\left|u_0^L(x_1,x_2)\right|\le C|2\pi L{|}^{-r-1/2},(x_1,x_2)\in {\Omega }_H^{2\pi }.$$

$$\int_{\Omega_{H}}\left[A_{p}\nabla u_{T}^{L}\cdot\nabla\overline{v}-k^{2}c_{p}u_{T}^{L}\overline{v}\right]\,\mathrm{d}x\,-\int_{\Gamma_{H}}T^{+}\left(u_{T}^{L}\big{|}_{\Gamma_{H}}\right)\overline{v}\,\mathrm{d}s\,=\int_{\Gamma_{H}}f_{L}\overline{v}\,\mathrm{d}s\,,$$

$$A_p(x)=\left|\det \nabla {\Phi }_p(x)\right|\left[{\left(\nabla {\Phi }_p(x)\right)}^{-1}{\left({\left(\nabla {\Phi }_p(x)\right)}^{-1}\right)}^{\top }\right]\in L^{\infty }(D,{\mathbb{R}}^{2\times 2});$$

$$c_p(x)=\left|\det \nabla {\Phi }_p(x)\right|\in L^{\infty }(D).$$

$$\int_{\Omega_{H}}\left[A_{p}\nabla u_{d}^{L}\cdot\nabla\overline{v}-k^{2}c_{p}u_{d}^{L}\overline{v}\right]-\int_{\Gamma_{H}}T^{+}\left(u^{L}_{d}\big{|}_{\Gamma_{H}}\right)\overline{v}\,\mathrm{d}s\,=\widetilde{b}(u_{0}^{L},v)$$

Full text

By transforming the original scattering problem into one defined in an infinite rectangle, the information of the surface is included in the coefficients. The numerical scheme contains two steps. The first step is to obtain an initial guess, i.e., the locations of both the periodic surfaces and the local perturbations, from a sampling method. The second step is to reconstruct the surface. As is proved in this paper, for some incident fields, the corresponding scattered fields carry little information of the perturbation. In this case, we use this scattered field to reconstruct the periodic surface. Then we could apply the data that carries more information of the perturbation to reconstruct the local perturbation. The Newton-CG method is applied to solve the associated optimization problems. Numerical examples are given at the end of this paper to show the efficiency of the numerical method.

1 Introduction

In this paper, we introduce the numerical method of the inverse scattering problem from a locally perturbed periodic surface. Both the periodic part and the local perturbation of the surface are unknown. The aim of the inveres problem is to reconstruct both of them from the near-filed measurement data.

Since the periodic surface is perturbed, the classical framework for the quasi-periodic scattering problems (i.e., quasi-periodic incident fields with periodic domains) no longer works. An efficient way to solve these problems is to apply the Floquet-Bloch transform. With the help of this Fourier-like transform, the original problem, which is defined in a 2D unbounded domain, is written into a new one defined in a 3D bounded domain. This method has been applied to perturbed periodic structures in [Coa12] and waveguide problems in [HN16]. For scattering problems with non-periodic incident fields and periodic surfaces, we refer to [LZ17a, LZ17c]. For problems with locally perturbed periodic surfaces, see [Lec17, LZ17b]. In the paper [Zha18], a high order numerical method has been proposed based on the Floquet-Bloch transform and this method is used in this paper to produce the measured data. For a fast imaging method to reconstruct the local perturbations in periodic media with the help of the Bloch transform, we refer to [CHN18].

The work in this paper is an extension to the joint work of the second author with Prof. Armin Lechleiter in [LZ18], where the periodic surface is assumed to be already known. A numerical method has been proposed to find out both the location and the shape of the local perturbation. The sampling method introduced by Ito, Jin and Zou (see [IJZ12]) was extended to find out the location, and a Newton-CG method was applied to reconstruct the shape. However, the setting in this paper is more difficult, i.e., the periodic surface is no longer known. Thus we have to find out the location without a known periodic surface, and also reconstruct both the periodic surface and the perturbation. In this case, the sampling method in [LZ18] does not work any more, which makes the problem much more challenging.

In this paper, we develop a numerical method for the inverse problem. The first task is to find out the locations of both the perturbation and the periodic surface. Since the previous sampling method does not work, we apply the rough surface reconstruction algorithm introduced in [LZZ18] to obtain the locations, when the perturbation is assumed to be existing in a relatively large domain. Then we apply the Newton’s method to reconstruct the shapes of both the periodic surface and the local perturbation. The reconstruction contains two steps. The first step is to reconstruct the periodic surface. From an estimate of the difference between field with and without perturbation, for certain incident fields, the scattered field with perturbation could be a good approximation of the one without perturbation. Thus in this case, the measured scattered field could be adopted to reconstruct the periodic surface. Base on the former approximation of the periodic surface, we apply the method in [LZ18] to find out the approximation of the local perturbation.

The rest of the paper is organized as follows. In Section 2, we recall the mathematical model of the direct scattering problem and the Floquet-Bloch transform based formulation. In Section 3, the estimation is considered for the difference between scattered fields with and without perturbation. The inverse problem is formulated in Section 4, and the Fréchet derivative and its adjoint operator are studied. In Section 5, we conclude the algorithm for inverse problems, including the sampling method for the initial guess and the iterative method for the reconstruction. In Section 6, we present two numerical results obtained from our algorithm.

2 Direct Scattering Problem

2.1 Mathematical Model

Given a bounded $2\pi$-periodic function $\zeta$, it defines a periodic surface

[TABLE]

Let the function $p$ be a compactly supported perturbation. For simplicity, suppose ${\rm supp}(p)\subset(-\pi,\pi)+2\pi J$, where $J\in\mathbb{Z}$ is an integer. Let $\zeta_{p}:=\zeta+p$ be the perturbed function and define

[TABLE]

Let the domain above $\Gamma$ be $\Omega$ and that above $\Gamma_{p}$ be $\Omega_{p}$.

Remark 1**.**

For simplicity, from Section 2 to Section 3, we fix $J=0$ for theoretical arguments.

In this paper, we assume that the surface $\Gamma_{p}$ is sound-soft. Given an incident field $u^{i}$ that satisfies $\Delta u^{i}+k^{2}u^{i}=0$ in $\Omega_{p}$, then it is scattered by $\Gamma_{p}$ and generates the scattered field $u^{s}$ (or equivalently, the total field $u=u^{i}+u^{s}$). For the mathematical model we refer to Figure 1. First, $u$ satisfies

[TABLE]

Second, as the surface $\Gamma_{p}$ is sound-soft,

[TABLE]

Moreover, the scattered field $u^{s}$ is propagating upwards. The Upward Propogation Radiation Condition (UPRC) is typically written as a double layer potential, see [CWZ98], and an alternative definition was introduced in [CM05, CE10]. Let $H$ be a real number that is larger than $\|\zeta\|_{\infty}$ and $\|\zeta_{p}\|_{\infty}$. Then the UPRC is written as

[TABLE]

where $\widehat{u}^{s}(\xi,H)$ is the Fourier transform of $u^{s}(x_{1},H)$. Define the Dirichlet-to-Neumann map $T^{+}$ by

[TABLE]

Let $\Gamma_{H}:=\left\{(x_{1},H):\,x_{1}\in\mathbb{R}\right\},$ then the UPRC is equivalent to

[TABLE]

Define the domain $\Omega^{p}_{H}:=\mathbb{R}\times(-\infty,H)\cap\Omega_{p}$ and we consider the problems (1)-(3) in the weighted Sobolev space $H_{r}^{1}(\Omega^{p}_{H})$, where the space $H_{r}^{1}(\Omega^{p}_{H})$ is defined by

[TABLE]

$\widetilde{H}_{r}^{1}(\Omega^{p}_{H})$ is the subspace of $H_{r}^{1}(\Omega^{p}_{H})$ such that all the elements vanish on $\Gamma_{p}$. Similarly, we can define the weighted spaces $H^{1/2}_{r}(\Gamma_{H})$ and $H^{-1/2}_{r}(\Gamma_{H})$. From [CE10], the operator $T^{+}$ is bounded and continuous from $H_{r}^{1/2}(\Gamma_{H})$ to $H_{r}^{-1/2}(\Gamma_{H})$ for all $|r|<1$.

The weak formulation of the scattering problem (1)-(3) is to find $u\in\widetilde{H}_{r}^{1}(\Omega^{p}_{H})$ such that

[TABLE]

for all $v\in\widetilde{H}^{1}_{r}(\Omega^{p}_{H})$ with compact support in $\overline{\Omega^{p}_{H}}$. The unique solvability of the variational problem (5) has been proved in [CE10]:.

Theorem 2**.**

For $|r|<1$, given any incident field $u^{i}$ and the function $f$ defined by (3) belongs to the space $H_{r}^{-1/2}(\Gamma_{H})$, the variational problems (5) has a unique solution $u\in\widetilde{H}^{1}_{r}(\Omega^{p}_{H})$.

Remark 3**.**

Although the unique solvability is proved for bounded surfaces, in this paper, the functions $\zeta$ and $\zeta_{p}$ are assumed to be at least Lipschitz continuous.

2.2 Floquet-Bloch transform

During the numerical process of the inverse problem, a Floquet-Bloch transform based numerical method is applied to solve the direct scattering problems. Thus in this section, we give a brief introduction to this method. Let $h$ and $H$ be two real numbers such that $h<\min\{\zeta,\zeta_{p}\}<\max\{\zeta,\zeta_{p}\}<H$, then define $D:=\mathbb{R}\times(h,H)$. Define the periodic cell $W$ and its dual-cell ${W^{*}}$ by

[TABLE]

Then let $D^{2\pi}=D\cap W\times\mathbb{R}$, $\Gamma^{2\pi}_{h}=W\times\{h\}$, $\Gamma_{H}^{2\pi}=W\times\{H\}$. Define the Bloch transform with period $2\pi$ in $D$ by

[TABLE]

Define the function space $H_{0}^{r}({W^{*}};H^{s}_{\alpha}(D^{2\pi}))$ by the closure of $C_{0}^{\infty}({W^{*}}\times D^{2\pi})$ with the following norm for $r\in\mathbb{N}$:

[TABLE]

The definition is extended to all $r\geq 0$ by interpolation between Hilbert space, and to $r\in\mathbb{R}$ by duality arguments. The property of the Bloch transform has been investigated in [LZ17b].

Theorem 4**.**

The Bloch transform is an isomorphism between $H_{r}^{s}(D)$ and $H_{0}^{r}({W^{*}};H^{s}_{\alpha}(D^{2\pi}))$. Further, when $s=r=0$, $\mathcal{J}_{D}$ is an isometry with the inverse

[TABLE]

and the inverse transform equals to the adjoint operator of $\mathcal{J}_{D}$.

Now we apply the Floquet-Bloch transform to the scattering problem (5). Following [LZ17b], the first task is to transform the original problem, which is defined in the non-periodic domain $\Omega_{H}^{p}$, to a periodic domain. In this paper, we choose $D$ as the periodic domain. Let $H_{0}$ be a real number that lies in the interval $(\min\{\zeta,\zeta_{p}\},H)$ and define the following two diffeomorphisms for $x\in\Omega^{p}_{H_{0}}$:

[TABLE]

Then extend them by the identity operator for $x_{2}\geq H_{0}$. From the assumption that ${\rm supp}(\zeta_{p}-\zeta)\subset W$, ${\rm supp}(\Phi_{\zeta}-\Phi_{\zeta_{p}})\subset D^{2\pi}$.

Let $u_{D}=u\circ\Phi_{\zeta_{p}}$, it is easily checked that $u_{D}$ satisfies the following variational equation:

[TABLE]

for all $v_{D}=v\circ\Phi_{\zeta_{p}}\in\widetilde{H}^{1}(D)$, where

[TABLE]

We define the matrix $A_{\zeta}$ and $c_{\zeta}$ by $\Phi_{\zeta}$ in the similar way, i.e.,

[TABLE]

As ${\rm supp}(\Phi_{\zeta}-\Phi_{\zeta_{p}})\subset D^{2\pi}$, the supports of both $A_{\zeta_{p}}-A_{\zeta}$ and $c_{\zeta_{p}}-c_{\zeta}$ are subsets of $D^{2\pi}$. Let $w:=\mathcal{J}_{D}u_{D}$, then it satisfies

[TABLE]

where

[TABLE]

Following the arguments in [Lec17, LZ17b], it is easy to prove that when the functions $\zeta$ and $\zeta_{p}$ are Lipschitz continuous, the variational problem (8) is equivalent to (5). When (5) has a unique solution of $\widetilde{H}_{r}^{1}(\Omega^{p}_{H})$ for some $|r|<1$, the problem (8) has a unique solution in $H_{0}^{r}({W^{*}};\widetilde{H}^{1}_{\alpha}(D^{2\pi}))$. Moreover, if the incident field $u^{i}\in H_{r}^{2}(\Omega^{p}_{H})$ and the surfaces are $C^{2,1}$, then the solution belongs to the space $H_{0}^{r}({W^{*}};\widetilde{H}^{2}_{\alpha}(D^{2\pi}))$. In [LZ17b], a convergent numerical method based on (8) has been proposed for the numerical solution, and a high order method has been proposed in [Zha18].

Remark 5**.**

The information of the periodic function $\zeta$ is included in $A_{\zeta}$ and $c_{\zeta}$, and the information of $p$ is included in $A_{\zeta_{p}}$ and $c_{\zeta_{p}}$. During the iteration process, when $\zeta$ and $p$ are updated, the matrices $A_{\zeta},\,A_{p}$ and the functions $c_{\zeta},\,c_{p}$ are updated. Thus we do not need to change the meshes during this process.

3 Approximation of the scattering problems with periodic surfaces

This section considers the difference between the scattered fields with and without the local perturbation. Let $u_{0}$ be the total field with the same incident field $u^{i}$ and periodic surface $\Gamma$, then $u_{0}$ satisfies the variational equation:

[TABLE]

From Theorem 2, if $f\in H_{r}^{-1/2}(\Gamma_{H})$ for some $|r|<1$, then the solution $u_{0}\in\widetilde{H}_{r}^{1}(\Omega_{H})$. In this paper, we assume that $r\in(0,1)$. As $u_{0}\in H_{r}^{1}(\Omega_{H})$, there is a constant $C$ that does not depend on $u_{0}$ such that

[TABLE]

We apply the translation to the first variable, i.e., to replace $x_{1}$ by $x_{1}+2\pi L$ for some $L\in\mathbb{Z}\setminus\{0\}$, and let $u^{i}_{L}(x_{1},x_{2}):=u^{i}(x_{1}+2\pi L,x_{2})$ be the incident field. As the surface is $2\pi$-periodic, the total field with the incident field $u_{L}^{i}$, denoted by $u_{0}^{L}$, is actually the function $u_{0}(x_{1}+2\pi L,x_{2})$. $u_{0}^{L}$ satisfies the following variational equation

[TABLE]

with $f_{L}(x_{1},x_{2}):=f(x_{1}+2\pi L,x_{2})$ on $\Gamma_{H}$. As

[TABLE]

the following estimate holds:

[TABLE]

Let $u^{L}$ be the solution of (5) with $f$ be replaced by $f_{L}$. Similar to the previous section, we can define a diffeomorphism $\Phi_{p}$ that maps $\Omega^{p}_{H}$ to $\Omega_{H}$ and $\Phi_{p}-I_{2}$ is supported in $\Omega_{H}\cap W\times\mathbb{R}$. Let $u_{T}^{L}:=u^{L}\circ\Phi_{p}$, it is easily checked that $u^{L}_{T}$ satisfies

[TABLE]

where

[TABLE]

Moreover, ${\rm supp}(A_{p}-I_{2}),\,{\rm supp}(c_{p}-1)\subset\Omega_{H}\cap W\times\mathbb{R}(:=\Omega^{2\pi}_{H})$. Then the difference $u_{d}^{L}:=u_{T}^{L}-u_{0}^{L}$ satisfies the following variational equation, i.e.,

[TABLE]

for any $v\in H^{1}(\Omega_{H})$ with compact support, where

[TABLE]

From the representation of $\widetilde{b}(\cdot,\cdot)$, it is a bounded sesquilinear form satisfies

[TABLE]

where $C$ is the constant depends only on $\zeta$ and $\zeta_{p}$. Thus the right hand side of (10) satisfies

[TABLE]

From the equivalence between (5) and (7), the equation (10) is uniquely solvable in $H^{1}(\Omega_{H})$ when the right hand side is a antilinear functional on $H^{1}(\Omega_{H})$. Thus

[TABLE]

Based on the above analysis, the total field $u_{0}^{L}$ is a good approximation of $u^{L}_{T}$ if $L$ is sufficiently large. Especially, let $\sigma$ be the noise level of the measured data. When $L\in\mathbb{Z}$ has a large enough absolute value such that $C|2\pi L|^{-r-1/2}<\delta$, $u^{L}_{T}$ could be treated as the “exact solution” of the non-perturbed periodic surface with the incident field $u_{T}^{L}$. Let

[TABLE]

then $\widetilde{u}^{L}_{T}$ is a good approximation of $u_{0}$. In this case, the solution $\widetilde{u}^{L}_{T}$ could be applied in the inverse problems to reconstruct the periodic surface.

4 Inverse Problem and the Newton-CG Method

The aim of the inverse problem is to reconstruct the unknown function $\zeta_{p}$ from the measured scattered data. The measured scattered field $U$ on $\Gamma_{H}$ is defined as

[TABLE]

where $\sigma$ is some noise added to the scattered data.

In the following, we always assume that $\zeta\in C^{2,1}(\mathbb{R})$ is a $2\pi$-periodic function and $p\in C^{2,1}(\mathbb{R})$ is a function that is compactly supported in $W+2\pi J$ for some $J\in\mathbb{Z}$.

Remark 6**.**

For the inverse problem, $J$ is an unknown integer and one task for the inverse problem is to find out the exact value of $J$. As is explained later, the integer $J$ could be found out by a sampling method (see [LZZ18]). Thus in this section, we treat it as a known one. For any $J\neq 0$, we can simply apply the translation $x\mapsto x-2\pi J$ to move the perturbation to the center of the domain (i.e., $J=0$). Thus for simplicity, we still assume that $J=0$ in this section.

Define the spaces

[TABLE]

In the following, we assume that $(\zeta,p)\in X\times Y$ and $\zeta_{p}:=\zeta+p$. The inverse problem is to find out $(\zeta,p)\in X\times Y$ such that the scattered field corresponding to $(\zeta,p)$ is the best approximation of $U$.

4.1 Scattering operator and its properties with respect to rough surfaces

We recall the inverse scattering problems from rough surfaces introduced in [CWP02]. Let $BC^{1,1}(\mathbb{R})$ be the space of bounded, Lipschitz continuous function.

Remark 7**.**

It is easily checked that $X,\,Y\subset BC^{1,1}(\mathbb{R})$.

Suppose $f\in BC^{1,1}(\mathbb{R})$ and the surface $\Gamma_{f}$ is defined by $f$. We can also define the domain $\Omega_{f}$ by the domain above $\Gamma_{f}$, and $\Omega^{f}_{H}$ by the domain between $\Gamma_{f}$ and $\Gamma_{H}$, where $H$ is a real number that is larger than $\|f\|_{\infty}$. Given an incident field $u^{i}$, we define the following scattering operator

[TABLE]

Then the inverse problem can be written as the optimization problem, i.e., to find $f\in BC^{1,1}(\mathbb{R})$ such that

[TABLE]

Let

[TABLE]

then the inverse problem is to find out the minimizer of the functional $F$ in the domain $BC^{1,1}(\mathbb{R})$. To solve the minimization problem, we have to study the properties of the scattering operator $S$ first.

Theorem 8**.**

The operator $S$ is differentialble, and its derivative $DS$ is represented as

[TABLE]

where $u^{\prime}\in H_{r}^{1}(\Omega^{f}_{H})$ satisfies

[TABLE]

Here $u$ is the total field of the scattering problem (1)-(3).

For the proof of this theorem we refer to [Kir93, CWP02].

For the Newton’s method, we also need the adjoint operator of the Fréchet derivative $DS$, which is explained in the following Theorem.

Theorem 9**.**

The adjoint operator of $DS(f)$, denoted by $\left[DS(f)\right]^{*}$ is given by

[TABLE]

where $\nu$ is the normal derivative upwards, $u$ is the total field and $z$ satisfies

[TABLE]

Remark 10**.**

During the iteration steps, the problems (17)-(19) and (21)-(23) will be solved several times. We can always apply the method introduced in Section 2.2 to transform the problems first into the one defined in the unbounded rectangle $D$ by the transform $\Phi_{\zeta_{p}}$, and then apply the Floquet-Bloch transform to obtain the new problem defined in the bounded domain ${W^{*}}\times D^{2\pi}$. For details of the solution of (21)-(23) we refer to Remark 12 in [LZ18].

4.2 Discretization for locally perturbed periodic surfaces

Let $\{\varphi_{1},\varphi_{2},\dots,\varphi_{M},\dots\}$ be a basis in the space $X$ and $\{\psi_{1},\psi_{2},\dots,\psi_{N},\dots\}$ be a basis in the space $Y$. For the positive integers $M$ and $N$, define the finite-dimensional subspaces of $X$ and $Y$ by:

[TABLE]

For the coefficients ${\bm{C}^{M}}=(c_{1}^{M},\dots,c_{M}^{M})\in\mathbb{R}^{M}$ and ${\bm{D}^{N}}=(d_{1}^{N},\dots,d_{N}^{N})\in\mathbb{R}^{N}$, then the elements $\zeta_{M}\in X_{M}$ and $p_{N}\in Y_{N}$ could be written as

[TABLE]

For the function $\zeta\in X$, there is a ${\bm{C}^{M}}\in\mathbb{R}^{M}$ such that $\zeta_{M}$ is the approximation in $X_{M}$ of $\zeta$. The argument also holds for $p_{N}\in Y_{N}$ and $p\in Y$.

Define the operators $A$ and $B$ by

[TABLE]

Then define the operator

[TABLE]

which maps the coefficients of both the periodic function and the local perturbation to the scattered field. Then we can define the functional $F$ in the finite dimensional space $\mathbb{R}^{M}\times\mathbb{R}^{N}$ by

[TABLE]

and the inverse problem is formulated by the following finite dimensional problem:

**Discrete Inverse Problem: ** to find $({\bm{C}^{M}},{\bm{D}^{N}})\in\mathbb{R}^{M}\times\mathbb{R}^{N}$ such that

[TABLE]

We apply the Newton-CG method to solve the descritized inverse problem. The linearized equation is

[TABLE]

where $\delta{\bm{C}^{M}}=(\delta c_{1}^{M},\dots,\delta c_{M}^{M})\in\mathbb{R}^{M}$ and $\delta{\bm{D}^{N}}=(\delta d_{1}^{N},\dots,\delta d_{N}^{N})\in\mathbb{R}^{N}$, $(DP)({\bm{C}^{M}},{\bm{D}^{N}})$ is the Fréchet derivative of $P$ at $({\bm{C}^{M}},{\bm{D}^{N}})$. Define

[TABLE]

then the linearized equation is written as

[TABLE]

First, we have to calculate the derivative of $P$. As an operator defined in the finite dimensional space $\mathbb{R}^{M}\times\mathbb{R}^{N}$, from direct calculation,

[TABLE]

Thus

[TABLE]

Given any $\delta{\bm{C}^{M}}\in\mathbb{R}^{M}$ and $\varphi\in L^{2}(\Gamma_{H})$,

[TABLE]

Let $Q=\left[(DS)(A({\bm{C}^{M}})+B({\bm{D}^{N}}))\right]^{*}$, then

[TABLE]

Similarly, we can also get

[TABLE]

In the numerical implementation, we solve the discrete inverse problem separately, i.e., first fix ${\bm{D}^{N}}$ and solve the minimization problem (25) to find out the solution ${\bm{C}^{M}}$.Then we fix ${\bm{C}^{M}}$ and solve the problem with respect to ${\bm{D}^{N}}$.

To solve the minimization problems we apply the Newton-CG method. To minimize the function $F({\bm{C}^{M}},{\bm{D}^{N}})$ with fixed ${\bm{D}^{N}}$, we apply the following Newton-CG method.

Similarly, we can also minimize the function $F({\bm{C}^{M}},{\bm{D}^{N}})$ with fixed ${\bm{C}^{M}}$ by the following algorithm:

5 Numerical implementation

5.1 Sampling method

In this section, we use the sampling method introduced in [LZZ18] to give an initial guess of the perturbed periodic surface, especially for the first term ${\bm{c}_{1}^{0}}$ of ${\bm{C}^{0}}$ and the integer $J$ of the perturbation.

Suppose that location $y$ of incident point sources $u^{i}(x,y)$ is on a horizontal line $\Gamma_{H}:=\{(x_{1},H)\leavevmode\nobreak\ |\leavevmode\nobreak\ x_{1}\in\mathbb{R}\}$ above the surface, we measure scattered Cauchy data $(u^{s},\partial_{\nu}u^{s})$ generated by these point sources and the perturbed periodic surface on $\Gamma_{H}$. Here, $\partial_{\nu}u^{s}$ denotes the normal derivative of $u^{s}$ on $\Gamma_{H}$ with the direction $(0,1)$.

We introduce the following imaging function

[TABLE]

where $y^{\prime}=(y_{1},-y_{2})$ and $z^{\prime}=(z_{1},-z_{2})$. From the analysis in [LZZ18], we can expect that the imaging function $I(z)$ takes a large value when $z\in\Gamma_{p}$ and decays as $z$ moves away from $\Gamma_{p}$. In this way, we give an initial guess of the perturbed surface.

In numerical computation, we choose $2P+1$ incident point sources which are located at $y_{j}=(jh_{inc},H),$ $j=-P,...,0,...P,$ here $h_{inc}$ is a fixed interval between two adjacent points. The measurement line $\Gamma_{H}$ is truncated to be ${\Gamma_{H,A}}:=\{x\in\Gamma_{H}\leavevmode\nobreak\ |\leavevmode\nobreak\ |x_{1}|<A\}$ which will be discretized uniformly into $2Q$ subintervals so the step size is $h_{mea}=A/Q$. In addition, the lower-half circle $\mathbb{S}_{-}$ in the second integral in (30) will also be uniformly discretized into $R$ grids with the step size $\Delta\theta=\pi/R$. Then for each sampling point $z$ we get the following discrete form of (30)

[TABLE]

Here, the measurement points are denoted by $x_{i}=(-A+ih,H),$ $i=0,1,...,2Q,$ and the normal directions are denoted by $d_{k}=(\sin(-\pi+k\Delta\theta),\cos(-\pi+k\Delta\theta)),$ $k=0,1,...,R.$

Suppose the sampling area is a rectangle denoted by $[a,b]\times[c,d]$. We set the numbers of sampling points in $x_{1}$-direction and $x_{2}$-direction to be $M_{1}$ and $M_{2}$, respectively. Then by (31), we get the indicator matrix $\left\{I_{A}(z_{ij})\right\}_{M_{1}\times M_{2}}$. For each $j$th-row of this matrix, we figure out the element with the largest value $I_{A}$ and denote the corresponding index by ${{max}}_{j}$. The initial guess for the first term ${\bm{c}^{0}_{1}}$ of $C^{1}$ can be deduced by the following formula

[TABLE]

5.2 Iteration method

From the last subsection, we have already decided the integer $J$. By translation on the first variable, i.e., to let $x_{1}$ be replaced with $x_{1}-2\pi J$, the perturbation is moved to $W$ (i.e., $J=0$), then the method in Section 4 could be applied to the reconstruction.

From Section 3, for an incident field $u^{i}\in H_{r}^{1}(\Omega^{p}_{H})$ for some $r\in(0,1)$, the measured data $U_{L}$ with the incident field $u^{i}(\cdot+2\pi L,\cdot)$ for $L\in\mathbb{Z}\setminus\{0\}$, could be applied to reconstruct the periodic function $\zeta$. The measured data with incident field $u^{i}$, denoted by $U_{0}$, is then applied to reconstruct the local perturbation $p$. So we conclude the algorithm for the inverse scattering problem.

6 Numerical results

In this section, we present two examples for the numerical method. We define two different periodic surfaces and local perturbations:

[TABLE]

We apply Algorithm 3 to the following two examples (see Figure 2):

Example 1. The periodic surface $\Gamma$ is defined by $\zeta_{1}$ and the local perturbation is defined by $p_{1}$;

Example 2. The periodic surface $\Gamma$ is defined by $\zeta_{2}$ and the local perturbation is defined by $p_{2}$.

For both the incident point sources and Herglotz wave functions, the scattered data are collected on $\Gamma_{A,H}$ with $A=25\pi,H=3$ and it is divided into $2Q=1500$ subintervals with the step length $h_{mea}=\pi/300$. Let $u_{s}$ be the scattered data (either the scattered field or its normal derivative) on $\Gamma_{A,H}$, then the measured data is defined as:

[TABLE]

where $\sigma=5\%$ is the noise level and randn presents random numbers from the standard normal distribution.

6.1 Sampling method

For the sampling method, we choose the sampling area to be a rectangle as $[-20\pi,20\pi]\times[1.2,1.9]$. The number of sampling points in $x_{1}$-direction and $x_{2}$-direction are set to be $M_{1}=1600$ and $M_{2}=400$, respectively. For the first surface, we put 41 incident point sources at $y_{j}=(j\pi,3)$ with $j=-20,-19,\dots,20.$ For the second surface, we put 21 incident point sources at $y_{j}=(2j\pi,3)$ with $j=-10,-9,\dots,10$. The wavenumber is chosen to be $k=3$ for both examples.

Use the indicator function introduced in (31), we can get a rough reconstruction of the original perturbed periodic surfaces in Figure 3 and 4. Note that the red dash lines in (c) are boundaries of the periodic cells. In each figure, we first present the profile of the original surface. Then the reconstructed result is given directly by the indicator function $I_{A}(z)$. Finally, in order to give the initail guess of ${\bm{c}^{0}_{1}}$ and the integer $J$ of the perturbation, we try to find out the points $z_{max}$ which get the largest value $I_{A}(z)$ in each vertical line and plot them in the last position of each figure.

By the end of the sampling step, we give out the value $J$ and ${\bm{c}^{0}_{1}}$. Roughly speaking, $J$ represents the location of the perturbation while ${\bm{c}^{0}_{1}}$ gives the vertical location of the periodic surface. From Figure 3 and 4, the locations of the perturbations are easily obtained, i.e., $J=-3$ for Example 1 and $J=2$ for Example 2. The initial guess of ${\bm{c}^{0}_{1}}$ is computed due to (32). By staightward calcultaions, we get ${\bm{c}^{0}_{1}}=1.4987$ and ${\bm{c}^{0}_{1}}=1.5216$. Both of these two results are very good approximations of the constant terms of both $\zeta_{1}$ and $\zeta_{2}$.

6.2 Newton’s method

For the Newton’s method, the Herglotz wave function is applied as the incident field (see Figure 5), i.e.,

[TABLE]

where

[TABLE]

Remark 11**.**

We could not use the point source as the incident fields since the fundamental solution $\Phi(x,y)=\frac{\mathrm{i}}{4}H_{0}^{(1)}(k|x-y|)$ belongs to the space $H_{r}^{1}(\Omega^{p}_{H})$ only if $r<0$. From [LZ17b], the direct solver introduced in Section 2.2 does not converge.

The incident field $u^{i}\in H_{r}^{1}(\Omega^{p}_{H})$ for any $r\in(0,1)$. Let $L=4$, then we use two incident fields $u^{i}$ and $u^{i}_{L}:=u^{i}(\cdot+2\pi L,\cdot)$. Let $u^{s}$ and $u^{s}_{L}$ be the scattered fields corresponding to the incident fields $u^{i}$ and $u^{i}_{L}$, and $u$, $u_{L}$ be the corresponding total fields. From the estimation in Section 3, the error between $u_{T}^{L}:=u\circ\Phi_{p}$ and $u_{0}^{L}$, which is the total field with incident field $u^{i}_{L}$ and the periodic surface, is bounded by:

[TABLE]

Note that the noise level is $\sigma=5\%$, $u_{T}^{L}$ could be treated as a good approximation of $u_{0}^{L}$ when the constant $C$ is assumed to be not too large.

Then we apply Algorithm 3 to reconstruct the perturbation with the known values $J$ and ${\bm{c}_{0}^{1}}$ from the sampling method. The reconstructs for Example 1 and Example 2 are shown in Figure 6 and 7, respectively. From the left pictures of the two figures, the periodic surfaces are well reconstructed; based on the results for the periodic surfaces, we can also reconstruct the local perturbations very well.

Acknowlegdments

This paper is devoted to Professor Armin Lechleiter. We will never forget him as a talented mathematician, a supportive colleague, and a dear friend.