Inverse scattering reconstruction of a three dimensional sound-soft axis-symmetric impenetrable object^*

To be able to use the uniqueness result from theorem 3.1, one must first find the axis of symmetry of the obstacle. Our algorithm to find the axis of symmetry is based on theorems 3.2 and 3.3.

Theorem 3.2. If the axis of symmetry of Ω is the z-axis, then for any ϕ ∈ [0, π], both the real and imaginary parts of the far field u^∞(θ, ϕ) of Ω by the incident wave u^inc = e^{ik d⋅x} with d = (d₁, d₂, d₃) ≠ (0, 0, ±1) are symmetric with respect to θ = θ₀ and θ = 2π − θ₀, where θ₀ satisfies $\mathrm{cos}\left({\theta }_{0}\right)={d}_{1}/\sqrt{{d}_{1}^{2}+{d}_{2}^{2}}$, and $\mathrm{sin}\left({\theta }_{0}\right)={d}_{2}/\sqrt{{d}_{1}^{2}+{d}_{2}^{2}}$. Here θ ∈ [0, 2π) is the azimuthal angle in the xy-plane from the positive x-axis and ϕ ∈ [0, π] is the altitude angle from the positive z-axis. If d = (0, 0, ±1), then the far field is axis-symmetric with respect to the z-axis.

Proof. For d ≠ (0, 0, ±1), without loss of generality, we may assume d₂ = 0, in which case we need to show that the far field is symmetric with respect to θ = 0 and θ = π. In fact, when d₂ = 0, we have ${u}^{\mathrm{i}\mathrm{n}\mathrm{c}}={\mathrm{e}}^{\mathrm{i}k\left({d}_{1}r\mathrm{cos}\theta +{d}_{3}z\right)}={\mathrm{e}}^{\mathrm{i}k\left({d}_{1}r\mathrm{cos}\left(-\theta \right)+{d}_{3}z\right)}$, so u^inc is symmetric with respect to θ = 0 and θ = π for any fixed z. On the other hand, the obstacle Ω is also symmetric with respect to the plane that is cut by θ = 0 and θ = π. By the uniqueness theorem of the exterior problem, we have that the scattered field u^scat is also symmetric with respect to the plane where θ = 0 or θ = π. The far field pattern can be obtained by using the equation

$$\begin{equation}{u}^{\infty }\left(\theta ,\phi \right)=\frac{1}{4\pi }{\int }_{\partial {\Omega}}\left\{{u}^{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}}\left(\mathbf{y}\right)\frac{\partial {\mathrm{e}}^{-\mathrm{i}k\hat{\mathbf{x}}\cdot \mathbf{y}}}{\partial \boldsymbol{\nu }\left(\mathbf{y}\right)}-\frac{\partial {u}^{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}}\left(\mathbf{y}\right)}{\partial \boldsymbol{\nu }\left(\mathbf{y}\right)}{\mathrm{e}}^{-\mathrm{i}k\hat{\mathbf{x}}\cdot \mathbf{y}}\right\}\mathrm{d}s\left(\mathbf{y}\right),\end{equation} \tag{ 3.1 }$$

where $\hat{\mathbf{x}}=\left(\mathrm{cos}\enspace \theta \enspace \mathrm{sin}\enspace \phi ,\mathrm{sin}\enspace \theta \enspace \mathrm{sin}\enspace \phi ,\mathrm{cos}\enspace \phi \right)$. Since all the components on the right-hand side of equation (3.1) are symmetric with respect to θ = 0 or θ = π, the far field must be symmetric with respect to θ = 0 or θ = π, too.

Similarly, when d = (0, 0, ±1), both the incident wave and the obstacle are axis-symmetric with respect to the z-axis, so is the scattered field and the far field pattern. □

Theorem 3.2 implies that the far field pattern is symmetric with respect to the axis of symmetry when the axis passes through the origin. For an axis-symmetric obstacle that is not centered at the origin, we have the following translation property for the far field.

Theorem 3.3. Let u^∞(θ, ϕ) be the far field pattern of Ω by the incident field u^inc(x) = e^{ik d⋅x} . For a shifted domain Ω_h := {x + h, x ∈ Ω} with a constant vector $h\in {\mathbb{R}}^{3}$, the far field ${u}_{h}^{\infty }\left(\theta ,\phi \right)$ generated by the incident wave u^inc(x) becomes

$$\begin{equation*}{u}_{h}^{\infty }\left(\theta ,\phi \right)={\mathrm{e}}^{\mathrm{i}k\left(\mathbf{d}-\hat{\mathbf{x}}\right)\cdot h}{u}^{\infty }\left(\theta ,\phi \right)\end{equation*}$$

where $\hat{\mathbf{x}}=\left(\mathrm{cos}\enspace \theta \enspace \mathrm{sin}\enspace \phi ,\mathrm{sin}\enspace \theta \enspace \mathrm{sin}\enspace \phi ,\mathrm{cos}\enspace \phi \right)$.

The proof of this theorem can be found in [26]. Readers are referred to [15, 16] for similar conclusions in elastic wave scattering.

Theorem 3.3 implies that the shifted domain of Ω simply changes the phase of the far field but not the modulus. Therefore, for a given far field data u^∞(θ, ϕ) of an unknown axis-symmetric obstacle Ω, |u^∞(θ, ϕ)| is the same as the modulus of the far field of Ω₀, where Ω₀ denotes a shifted Ω with its axis centered at the origin. Thus by theorem 3.2, up to a rotation, |u^∞(θ, ϕ)| is symmetric with respect to θ = θ₀ and θ = 2π −θ₀ for a given θ₀ and any ϕ ∈ [0, π]. This rotation angle is exactly the orientation of the axis of symmetry of the unknown obstacle Ω. Once the axis is parallel to the z-axis, we can make use of the phase information to determine the x and y coordinates of the axis. In particular, by theorem 3.3, if we multiply the far field u^∞(θ, ϕ) by an appropriate factor ${\mathrm{e}}^{\mathrm{i}k\hat{x}\cdot h}$, the real and imaginary parts of the new far field will be symmetric with respect to θ = θ₀ and θ = 2π −θ₀ for a given θ₀ and any ϕ ∈ [0, π].

To summarize, we propose a three-step method to locate the axis of symmetry of the obstacle. In the first step, we evaluate the far field pattern based on the measured scattered field ${\mathbf{u}}_{k,\mathbf{d}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}$ on the sphere $\partial \mathcal{B}$. Next, we determine the orientation of the axis of symmetry. In the third step, we find the location of its center. A detailed description of the algorithm follows:

• (a)  
  Step 1 (Evaluate the far field): for a fixed wavenumber k and d, collect the measured scattered field ${\mathbf{u}}_{k,\mathbf{d}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}\left(\mathbf{x}\right)$ on $\partial \mathcal{B}$ due to the incident plane wave u^inc(x). From the measurements of the scattered field on $\partial \mathcal{B}$, one can obtain the scattered field u^scat(x) anywhere on $\partial \mathcal{B}$ by using interpolation. To find the corresponding far field u^∞(θ, ϕ), we solve the boundary integral equation

  $$\begin{equation*}\left(\frac{1}{2}\mathcal{I}+{\mathcal{D}}_{\partial \mathcal{B}}+ik{\mathcal{S}}_{\partial \mathcal{B}}\right)\eta \left(\mathbf{x}\right)={u}_{k,\mathbf{d}}^{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}}\left(\mathbf{x}\right),\quad \text{on}\enspace \partial \mathcal{B}.\end{equation*}$$

  where ${\mathcal{D}}_{\partial \mathcal{B}}$ and ${\mathcal{S}}_{\partial \mathcal{B}}$ are the single and double layer potentials defined on $\partial \mathcal{B}$. Once η(x) is found, the far field can be evaluated using the formula

  $$\begin{equation*}{u}^{\infty }\left(\theta ,\phi \right)=\frac{1}{4\pi }{\int }_{\partial \mathcal{B}}\left\{\frac{\partial {\mathrm{e}}^{-\mathrm{i}k\hat{\mathbf{x}}\cdot \mathbf{y}}}{\partial \boldsymbol{\nu }\left(\mathbf{y}\right)}+ik{\mathrm{e}}^{-\mathrm{i}k\hat{\mathbf{x}}\cdot \mathbf{y}}\right\}\eta \left(\mathbf{y}\right)\mathrm{d}s\left(\mathbf{y}\right).\end{equation*}$$
• (b)  
  Step 2 (Determine the orientation): suppose the far field u^∞(θ, ϕ) is given at 0 = θ₀ < θ₁ < ⋯ < θ_i < ⋯ < θ_n = 2π for θ and 0 = ϕ₀ < ϕ₁ < ⋯ < ϕ_j < ⋯ < ϕ_m = π for ϕ. Check if |u^∞(θ, ϕ)| is symmetric on the horizontal cross section by taking each (θ_i , ϕ_j ) as the north pole. If that is found within a certain accuracy, we take (θ_i , ϕ_j ) as the orientation of the axis of symmetry.
• (c)  
  Step 3 (Determine the location): assume the orientation of the axis of symmetry is parallel to the z-axis. In order to find the x and y coordinates of the axis, we send two incident waves ${u}_{1}^{\mathrm{i}\mathrm{n}\mathrm{c}}\left(\mathbf{x}\right)={\mathrm{e}}^{\mathrm{i}kx}$ and ${u}_{2}^{\mathrm{i}\mathrm{n}\mathrm{c}}\left(\mathbf{x}\right)={\mathrm{e}}^{\mathrm{i}ky}$, and evaluate their far field data ${u}_{1}^{\infty }\left(\theta ,\phi \right)$ and ${u}_{2}^{\infty }\left(\theta ,\phi \right)$ respectively. Suppose the obstacle is located in the area [x_min, x_max] × [y_min, y_max]. By theorem 3.3, we can determine h₁ such that the real and imaginary parts of the shifted far field ${u}_{2}^{\infty }\left(\theta ,\phi \right){\mathrm{e}}^{\mathrm{i}k\mathrm{cos}\theta \mathrm{sin}\phi {h}_{1}}$ are symmetric with respect to θ = π/2 and θ = 3π/2 for any ϕ ∈ [0, π] by searching a uniform grid in the interval [x_min, x_max]. Similarly, there exists h₂ ∈ [y_min, y_max] such that the real and imaginary parts of the shifted far field ${u}_{1}^{\infty }\left(\theta ,\phi \right){\mathrm{e}}^{\mathrm{i}k\mathrm{sin}\theta \mathrm{sin}\phi {h}_{2}}$ are symmetric with respect to θ = 0 and θ = π for any ϕ ∈ [0, π]. In the end, we take (h₁, h₂) as the x and y coordinates of the center of the obstacle.

Once the axis of symmetry of the obstacle is obtained, we can reconstruct the shape of the obstacle. Our goal is to find an approximation of the simple open curve $\gamma :\left[0,1\right]\to {\mathbb{R}}^{2}$ that generates the surface of the obstacle using the measurements ${\mathbf{u}}_{{k}_{j}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}$ for j = 1, ..., N_k . The idea is to apply the RLA to solve a sequence of single frequency inverse scattering problems.

3.2.1. Inverse scattering problem for a single frequency data

Using single frequency data, we can recast the inverse problem as the optimization problem

$$\begin{equation}\tilde {\gamma }=\mathrm{arg}\enspace \underset{\gamma }{\mathrm{min}}{\Vert}{\mathbf{u}}_{k}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}-{\mathcal{F}}_{k}\left(\gamma \right){\Vert},\end{equation} \tag{ 3.2 }$$

where $\tilde {\gamma }$ is an approximation of the exact curve γ that generates the obstacle ∂Ω.

The problem (3.2) is both nonlinear and ill-posed. To treat the nonlinearity, we apply the iterative damped Gauss–Newton method. First, an initial guess for the approximation of the generating curve γ of the boundary ∂Ω is chosen, say γ₀. Next, in each step of this method, given an approximation γ_j of the generating curve at the jth step, we update the curve to obtain

$$\begin{equation*}{\gamma }_{j+1}\left(t\right)={\gamma }_{j}\left(t\right)+\alpha \delta \gamma \left(t\right),\quad t\in \left[0,1\right],\end{equation*}$$

with α > 0 being a chosen constant and δγ = (δγ_r , δγ_z ). Note that the curve γ_j in each step of the method generate an axis-symmetric obstacles with boundary ∂Ω_j . To obtain the update δγ, we solve

$$\begin{equation}{\partial }_{\gamma }{\mathcal{F}}_{k}\left({\gamma }_{j}\right)\delta \gamma ={\mathbf{u}}_{k}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}-{\mathcal{F}}_{k}\left({\gamma }_{j}\right),\end{equation} \tag{ 3.3 }$$

where ${\partial }_{\gamma }{\mathcal{F}}_{k}\left({\gamma }_{j}\right)=\left[{\partial }_{\gamma }{\mathcal{F}}_{k,{\mathbf{d}}_{1}}\left({\gamma }_{j}\right);\cdots \enspace ;{\partial }_{\gamma }{\mathcal{F}}_{k,{\mathbf{d}}_{{N}_{d}}}\left({\gamma }_{j}\right)\right]$ and ${\partial }_{\gamma }{\mathcal{F}}_{k,{\mathbf{d}}_{i}}\left({\gamma }_{j}\right)$, i = 1, ..., N_d are, respectively, the Fréchet derivatives of ${\mathcal{F}}_{k}$ and ${\mathcal{F}}_{k,{\mathbf{d}}_{i}}$ with respect to ∂Ω evaluated at ∂Ω_j . As shown in [14] and references within, once we obtain the normal derivative of the total field ${u}_{k,{\mathbf{d}}_{i}}$ by solving problem (1.2) for the obstacle Ω_j with the incoming plane wave ${u}_{k,{\mathbf{d}}_{i}}^{\mathrm{i}\mathrm{n}\mathrm{c}}={\mathrm{e}}^{\mathrm{i}k\mathbf{x}\cdot {\mathbf{d}}_{i}}$, the value of $v\left(\mathbf{x}\right)=\left[{\partial }_{\gamma }{\mathcal{F}}_{k,{\mathbf{d}}_{i}}\left({\gamma }_{j}\right)\delta \gamma \right]\left(\mathbf{x}\right)$ is obtained by solving the Helmholtz equation

$$\begin{equation}\begin{cases}v\left(\mathbf{x}\right)+{k}^{2}v\left(\mathbf{x}\right)=0\quad \text{for}\enspace \mathbf{x}\in {\mathbb{R}}^{3}{\backslash}{\overline{{\Omega}}}_{j},\quad \hfill \\ v\left(\mathbf{x}\right)=-\left(\delta {\Gamma}\left(\mathbf{x}\right)\cdot \boldsymbol{\nu }\left(\mathbf{x}\right)\right)\left[\frac{\partial {u}_{k,{\mathbf{d}}_{i}}}{\partial \boldsymbol{\nu }}\right]\left(\mathbf{x}\right)\quad \text{for}\enspace \mathbf{x}\in \partial {{\Omega}}_{j},\quad \hfill \end{cases}\end{equation} \tag{ 3.4 }$$

where v also satisfies the Sommerfeld radiation condition and δΓ(x) = (δγ_r cos θ, δγ_r  sin θ, δγ_z ) is the update vector for x = (r(t), θ, z(t)) ∈ ∂Ω_j . Since it holds

$$\begin{equation*}\delta {\Gamma}\cdot \boldsymbol{\nu }=\left(\delta {\gamma }_{r}\enspace \mathrm{cos}\enspace \theta ,\delta {\gamma }_{r}\enspace \mathrm{sin}\enspace \theta ,\delta {\gamma }_{z}\right)\cdot \left({\nu }_{r}\enspace \mathrm{cos}\enspace \theta ,{\nu }_{r}\enspace \mathrm{sin}\enspace \theta ,{\nu }_{z}\right)=\delta \gamma \cdot \nu ,\end{equation*}$$

by combining (3.3) and (3.4), we are able to obtain the update δγ.

The iterations are repeated until a stopping criteria is reached. The stopping criteria can be the total number of iterations N_it, the residual achieving ${\Vert}{\mathbf{u}}_{k}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}-{F}_{k}\left({\gamma }_{j}\right){\Vert}{\leqslant}{{\epsilon}}_{r}$, with _r > 0, the difference of the curve evaluated in a set of points between consecutive steps being smaller than a certain value _s > 0, or others. A summary of the damped Gauss–Newton method is presented in algorithm 1.

Algorithm 1. Damped Gauss–Newton method.

1: Input: scattered field measurements ${\mathbf{u}}_{k}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}$, initial guess γ0, parameters α, Nit , r , and s .

2: Set j = 0, γ = γ0 and γ−1 = γ0(1 + 2||γ0||−1 s ).

3: while j < Nit and ${\Vert}{\mathbf{u}}_{k}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}-{\mathcal{F}}_{k}\left(\gamma \right){\Vert}{ >}{{\epsilon}}_{r}$ and ||γj − γj−1|| > s do

4:  Calculate ${\mathcal{F}}_{k}\left({\gamma }_{j}\right)$ and ${\partial }_{\gamma }{\mathcal{F}}_{k}\left({\gamma }_{j}\right)$.

5:  Solve ${\partial }_{\gamma }{\mathcal{F}}_{k}\left({\gamma }_{j}\right)\delta \gamma ={\mathbf{u}}_{k}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}-{\mathcal{F}}_{k}\left({\gamma }_{j}\right)$.

6:  γj+1 ← γj + αδγ

7:  j ← j + 1

8: end while

As mentioned, problem (3.2) is highly ill-posed. Various ways were proposed to deal with the ill-posedness of the inverse scattering problem, including Tykhonov regularization, truncated SVD [14], the use of a bandlimited representation of the domain [6], etc. In this work, we choose to search for a bandlimited representation of the generating curve γ. We represent the generating curve γ and the update δγ as

$$\begin{equation}\gamma \left(t\right)=p\left(t\right)\left(\mathrm{cos}\left(\pi \left(t-0.5\right)\right),\mathrm{sin}\left(\pi \left(t-0.5\right)\right)\right)\end{equation} \tag{ 3.5 }$$

and

$$\begin{equation}\delta \gamma \left(t\right)=h\left(t\right)\left(\mathrm{cos}\left(\pi \left(t-0.5\right)\right),\mathrm{sin}\left(\pi \left(t-0.5\right)\right)\right)\end{equation} \tag{ 3.6 }$$

where p(t) and h(t) are given by

$$\begin{equation}p\left(t\right)={p}_{0}^{\mathrm{c}}+\sum _{j=1}^{{N}_{p}}\left({p}_{j}^{\mathrm{c}}\mathrm{cos}\left(2\pi j\left(t-0.5\right)\right)+{p}_{j}^{\mathrm{s}}\mathrm{sin}\left(2\pi j\left(t-0.5\right)\right)\right)\end{equation} \tag{ 3.7 }$$

and

$$\begin{equation}h\left(t\right)={h}_{0}^{\mathrm{c}}+\sum _{j=1}^{{N}_{p}}\left({h}_{j}^{\mathrm{c}}\mathrm{cos}\left(2\pi j\left(t-0.5\right)\right)+{h}_{j}^{\mathrm{s}}\mathrm{sin}\left(2\pi j\left(t-0.5\right)\right)\right)\end{equation} \tag{ 3.8 }$$

with ${p}_{j}^{\mathrm{c}}$, ${p}_{j}^{\mathrm{s}}$, ${h}_{j}^{\mathrm{c}}$ and ${h}_{j}^{\mathrm{s}}$ being constant coefficients for the jth cosine and sine modes, j = 1, ..., N_p . From Heisenberg's uncertainty principle for waves, we have that sub-wavelength features of the scatterer are present in the evanescent modes of the signal and are not detectable in finite precision. Consequently, the main advantage of this representation is that if we choose the bandlimit parameter ${N}_{p}=\mathcal{O}\left(k\right)$, the system of equations on (3.3) becomes well-conditioned. The second advantage of choosing this representation comes from the easy and fast evaluation of the polynomials p(t) and h(t) by using non-uniform FFT [20, 30]. The main disadvantage of choosing this representation stems from the fact that this representation is ideal for star-shaped figures. If the obstacle that we are trying to recover is not star-shaped, this representation will probably not work in terms of providing a high resolution reconstruction. An alternative is to use a bandlimited curve smoother like the one presented in [4, 6].

3.2.2. Inverse scattering problem using multiple frequency data

On the one hand, due to Heisenberg's uncertainty principle, the amount of information about the shape of the scatterer that can be stably recovered from measurements of the scattered field at frequency k is proportional to $\mathcal{O}\left(k\right)$. This means that smaller features of the obstacle that have magnitude proportional to the sub-wavelength spectrum are extremely difficult to recover using finite precision. On the other hand, there are also inherent limitations to Newton-type methods for inverse scattering at a single frequency. When the incident field has larger wavelength, a low-resolution approximation of the inhomogeneity can be obtained using a simple initial guess. For problems with incident field with smaller wavelength, the initial guess for the iterative method must be close to the solution for the method to converge. This interplay between obtaining a low resolution reconstruction for large wavelengths using a simple initial guess, and the need to have a very good initial guess when using small wavelengths, together with the natural limitation on the amount of information that can be obtained using single frequency data led to the proposal of the RLA [1, 10, 11].

In the RLA, a sequence of increasingly complicated nonlinear optimization problems like (3.2) at successively higher frequencies are solved using a continuation path in frequency. At a given frequency k, one uses the damped Gauss–Newton method to solve the inverse problem (3.2) and obtain an approximate solution γ^(k) to the curve γ. This solution is used as the initial guess for the damped Gauss–Newton method to solve the nonlinear optimization problem with scattered field data at frequency k + δk, where δk > 0 is sufficiently small. A summary of the RLA is presented in algorithm 2.

Algorithm 2. RLA with damped Gauss–Newton method.

1: Input: scattered field measurements ${\mathbf{u}}_{{k}_{j}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}$, for j = 1, ..., Nk and ${k}_{1}{< }{k}_{2}{< }\dots {< }{k}_{{N}_{k}}$, initial

guess γ0, parameters α, Nit, r , and s .

2: Set γ(0) = γ0.

3: for j = 1, ..., Nk do

4:  Apply the damped Gauss–Newton method with initial γ(j−1), parameters α, Nit , r , and s

to obtain an approximation $\tilde {\gamma }$ of the curve as a solution.

5:  ${\gamma }^{\left(k\right)}{\leftarrow}\tilde {\gamma }$.

6: end for

We test part 1 of our framework by determining the orientation and location of the axis of symmetry of an unknown obstacle. In particular, we are trying to determine the axis of an oblique ellipsoid generated by rotating and shifting from a standard one, whose parameterization of generating curve is given by

$$\begin{equation*}\gamma \left(t\right)=\left(x\left(t\right),y\left(t\right)\right)=\left(\mathrm{cos}\left(t\pi -0.5\right),2\mathrm{sin}\left(t\pi -0.5\right)\right),t\in \left[0,1\right].\end{equation*}$$

The axis of the ellipsoid is oriented at $\theta =\frac{\pi }{4}$, $\phi =\frac{\pi }{3}$ with the center located at (2, 2, 0), as shown in figure 2(a). Using the incident wave u^inc(x) = e^{ik x⋅d} with k = 3 and d = (cos(π) sin(π/8), sin(π) sin(π/8), cos(π/8)), we measure the scattered field u^scat(x) on $\partial \mathcal{B}$ with $\mathcal{B}$ being radius 10 and centered at the origin. The far field pattern ${u}^{\infty }\left(\hat{\mathbf{x}}\right)$ is found according to step 1 in part 1 of our algorithm. The modulus of ${u}^{\infty }\left(\hat{\mathbf{x}}\right)$ is shown in figure 2(b). In particular, the symmetry cannot be seen directly. We therefore test different rotation angles as the north pole as stated in the step 2 of part 1 by choosing m = n = 100. By applying this test, we successfully find the orientation angle θ = 0.25π and ϕ = 0.33π, which is very close to the exact solution. The far field pattern after rotation is shown in figure 2(c). We also plot the cross section of the far field pattern at ϕ = π/4 before and after the rotation in figure 3(a). One clearly sees that the symmetry is recovered if the correct rotation is found.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Once the orientation is found, our next step is to determine the location of the axis on the xy−plane. We collect the far field data by sending two incident waves, respectively. One is e^ikx and the other is e^iky . Since the location is not at the center, the far field patterns are not symmetric anymore. However, by finding the corresponding phase function e^iα and e^iβ from step 3 of part 1, with α = kh₁ cos θ sin ϕ, β = kh₂ sin θ sin ϕ and (h₁, h₂) = (2, 2), we are able to recover the symmetry of the far field pattern, as shown in figures 3(b) and (c) for a cross section at ϕ = π/4.

In this example, let the objective functional be ${f}_{k}\left(\gamma \right)={\Vert}{\mathbf{u}}_{k}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}-{\mathcal{F}}_{k}\left(\gamma \right){\Vert}$. We consider the case where the obstacle is a sphere of radius 1, with generating curve $\gamma :\left[0,1\right]\to {\mathbb{R}}^{2}$ given by $\gamma \left(t\right)=\left(\right.\mathrm{cos}\left(\pi \left(t-0.5\right),\mathrm{sin}\left(\pi \left(t-0.5\right)\right)\right)$. Applying the damped Gauss–Newton method one tries to recover at each step a polynomial p(t) = p₀. The shape reconstruction problem of finding p₀ turns into the problem of finding the root of a single variable nonlinear equation. To guarantee the convergence of the damped Gauss–Newton method, the initial guess must be in the same local convexity set of f_k as the solution. This example aims to illustrate the interplay between the wavenumber of the incident wave and the local set of convexity near the solution.

We use for each experiment one incident plane wave with incident direction d = (cos(π/9), 0, sin(π/9)) and wavenumber k_m , such that k₁ = 1, k_m = 2.5(m − 1), m = 2, ..., 13. The scattered data is measured at receptors ${x}_{\theta ,\phi }=\rho \left(\mathrm{cos}\left({\phi }_{j}\right)\mathrm{sin}\left({\theta }_{l}\right),\mathrm{sin}\left({\phi }_{j}\right)\mathrm{sin}\left({\theta }_{l}\right),\mathrm{cos}\left({\theta }_{l}\right)\right)$, with ρ = 10, θ_l = lπ/11, ϕ_j = 2jπ/10 for l, j = 1, ..., 10. Since we want to show the relation between the wavenumber and the local convexity set of f_k , we do not add noise to the measurements.

We calculate ${f}_{{k}_{m}}\left({\tilde {\gamma }}^{\left(j\right)}\right)$ at the curves ${\tilde {\gamma }}^{\left(j\right)}\left(t\right)={p}_{0}^{\left(j\right)}\left(\mathrm{cos}\left(t\right),\mathrm{sin}\left(t\right)\right)$, where ${p}_{0}^{\left(j\right)}=0.01j$, for j = 1, ..., 1000. Using the values of the objective functional, we are able to identify the local set of convexity in which p₀ = 1 is located. We denote this set to be $\left[a,b\right]$, where 0 < a is the largest point smaller than 1 where f_k attains a local maximum, and b is the smallest point larger than 1 where the function attains a local maximum.

In figure 4(a), we present the value of ${f}_{k}\left({\tilde {\gamma }}^{\left(j\right)}\right)$ for k = 1, 5, 10, 20 and 30. In figure 4(b), we plot three lines: the line for p₀ = 1, the line with values of b and the line with values of a. The values of a and b are also available in table 3 with the wavenumber and its respective wavelength λ = 2π/k.

Zoom In Zoom Out Reset image size

As expected, with the increasing value of the wavenumber k, the size of the interval $\left[a,b\right]$ decreases. Also, as the wavenumber increases, f_k presents multiple local minima, which shows the nonlinearity of the inverse problem.

In this example, we use multifrequency scattered data generated by four different geometric configurations of incident waves and receptors to recover the shape of an obstacle generated by the rotation around the z-axis of the curve $\gamma :\left[0,1\right]\to {\mathbb{R}}^{2}$, given by $\gamma \left(t\right)=p\left(t\right)\left(\mathrm{cos}\left(\pi \left(t-0.5\right)\right),\mathrm{sin}\left(\pi \left(t-0.5\right)\right)\right)$, where

$$\begin{equation*}p\left(t\right)=1.5+\left(\right.0.3\enspace \mathrm{cos}\left(8\pi \left(t-0.5\right)/0.5\right).\end{equation*}$$

A three-dimensional rendering of the obstacle can be seen in figure 5(a).

Zoom In Zoom Out Reset image size

We set up the incident wave directions and the receptors positions in four different geometric configurations:

• (a)  
  The incoming direction of the incident waves is d = (−1, 0, 0). The scattered field is measured at N_r = 100 receptors located at the points ${\mathbf{x}}_{m}=10\left(\mathrm{cos}\left({\theta }_{m}\right),0,\mathrm{sin}\left({\theta }_{m}\right)\right)$, with θ_m = −π/2 + mπ/(N_r − 1), m = 0, ..., N_r − 1. See figure 5(b).
• (b)  
  The incoming direction of the incident waves is d = (−1, 0, 0). The scattered field is measured at N_r = 100 receptors located at the points ${\mathbf{x}}_{m}=10\left(\mathrm{cos}\left({\theta }_{m}\right),0,\mathrm{sin}\left({\theta }_{m}\right)\right)$, with θ_m = −π/2 + ((N_θ − N_r )/2 + m)π/(N_θ − 1), m = 0, ..., N_r − 1, and N_θ = 3100. See figure 5(c).
• (c)  
  The incoming direction of the incident waves is d = (0, 0, −1). The scattered field is measured at N_r = 100 receptors located at the points ${\mathbf{x}}_{m}=10\left(\mathrm{cos}\left({\theta }_{m}\right),0,\mathrm{sin}\left({\theta }_{m}\right)\right)$, with θ_m = mπ/(N_r − 1), m = 0, ..., N_r − 1. See figure 5(b).
• (d)  
  The incoming direction of the incident waves is d = (0, 0, −1). The scattered field is measured at N_r = 100 receptors located at the points ${\mathbf{x}}_{m}=10\left(\mathrm{cos}\left({\theta }_{m}\right),0,\mathrm{sin}\left({\theta }_{m}\right)\right)$, with θ_m = ((N_θ − N_r )/2 + m)π/(N_θ − 1), m = 0, ..., N_r − 1 and N_θ = 3100. See figure 5(c).

In the configuration 1, the incident wave illuminates the obstacle in the direction perpendicular to the axis of symmetry and the receptors are located such that it is possible to see the entire obstacle from their positions (taking symmetry into consideration). In configurations 2, 3 and 4, the position of the receptors provide only limited information about the obstacle.

To obtain measurements, we start by computing the scattered field data ${u}_{k,\mathbf{d}}^{\text{scat}}$ at frequencies k_q = 0.5 + (q − 1)0.25, with q = 1, ..., 25. Next, to avoid inverse crimes, we add 2% noise to the measured scattered data, using the formula

$$\begin{equation}{\left({\mathbf{u}}_{k,\mathbf{d}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}\right)}_{m+1}={u}^{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}}\left({\mathbf{x}}_{m}\right)+0.02\frac{{\epsilon}}{{\Vert}{\epsilon}{\Vert}}\vert {u}_{k,\mathbf{d}}^{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}}\left({\mathbf{x}}_{m}\right)\vert \end{equation} \tag{ 4.1 }$$

where ${\left({\mathbf{u}}_{k,\mathbf{d}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}\right)}_{m+1}$ is the (m + 1)^th coordinate of the vector ${\mathbf{u}}_{k,\mathbf{d}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}$, m = 0, ..., N_r − 1, = ₁ + i₂, and ₁ and ₂ are chosen from the random normal distribution $\mathcal{N}\left(0,1\right)$ with mean zero and variance one.

We apply the RLA with the damped Gauss–Newton method at each frequency. We set the stopping criteria of the Gauss–Newton method to be the maximum number of iterations N_it = 10, the update step size should be no smaller than _s = 0.03 and the residual smaller than _r = 0.03. We also include as a stopping criteria any residual increase from one step to another. We set the damping parameter α = 0.1 for the damped Gauss–Newton method at frequency k = 0.5, and α = 0.1/||h|| for all other frequencies, where ||h|| is the 2-norm of the vector of coefficients of the update. We set the number of modes in the polynomial representing the update h to be N_p = ⌊2k⌋.

In figures 6(a)–(d), we present the reconstructions at frequency k = 6.5 for the configurations 1, 2, 3 and 4, respectively. In addition, in figure 7(a) 1, 7(b) 2, 7(c) 3 and 7(d), the cross section of the original obstacle and the reconstructions at k = 3.25 and k = 6.5 are shown for configurations 1, 2, 3 and 4, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The reconstruction obtained using configuration 1 is more accurate than the reconstructions obtained using the other configurations. This behavior was expected since the placement of the receptors in configuration 1 allows for obtaining information from a larger part of the obstacle.

This example is a continuation of example 3. We use the same scattered data with 2% noise that was generated for the configuration 1 in example 3 to recover the obstacle in figure 5(a). In example 3, at each frequency, the number of modes used to approximate the boundary of the obstacle and the update obtained by the Gauss–Newton step is N_p = ⌊2k⌋. Since the frequency varied from k = 0.5 to 6.5, the number of modes used varied from 1 to 13.

In example 4, instead of letting the number of modes increase freely with the frequency, we set it to be ${N}_{p}=\mathrm{min}\left\{{N}_{\mathrm{max}},\lfloor 2k\rfloor \right\}$, where we choose N_max = 8, 10 and 13. All the other parameters for both the RLA and the Gauss–Newton method are the same as in example 3.

In figures 8(a)–(c), we present the reconstructions at frequency k = 6.5 using 8, 10 and 13 modes, respectively. Figure 8(d) has the cross section of the original obstacle and the reconstructions using 8, 10 and 13 modes.

Zoom In Zoom Out Reset image size

As expected the reconstruction using 8 modes is more precise than the other reconstructions since the original obstacle only includes 8 modes. When we increase the number of modes used for the reconstruction, the results become increasingly worse due to the oscillations introduced by the higher order modes and the ill-posedness of the inverse problem. To deal with the oscillations introduced in the higher order modes, we introduce an appropriate filter in the next example.

In this example, we use multifrequency scattered data to reconstruct an obstacle with the shape of a land mine, see figure 9. One must use a very large number of modes N_p to recover the sharp corners of the obstacle based on a trigonometric representation.

Zoom In Zoom Out Reset image size

To generate the simulated scattered data ${\mathbf{u}}_{{k}_{j},{\mathbf{d}}_{\mathrm{s}}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}}$, we used incident plane waves given by ${u}_{{k}_{j},{\mathbf{d}}_{\mathrm{s}}}^{\mathrm{i}\mathrm{n}\mathrm{c}}\left(\mathbf{x}\right)={\mathrm{e}}^{\mathrm{i}{k}_{j}{\mathbf{d}}_{\mathrm{s}}\cdot \mathbf{x}}$, where k_j = 0.5 + (j − 1)0.25, j = 1, ..., 78, and

$$\begin{equation*}{\mathbf{d}}_{s}=\left(\mathrm{sin}\left(m\pi /6\right)\mathrm{cos}\left(2n\pi /5\right),\mathrm{sin}\left(m\pi /6\right)\mathrm{sin}\left(2n\pi /5\right),\mathrm{cos}\left(m\pi /6\right)\right),\end{equation*}$$

with s = (m − 1)5 + n, m, n = 1, ..., 5 and beyond that ${\mathbf{d}}_{26}=\left(0,0-1\right)$, and d₂₇ = (0, 0, 1). The scattered field is measured at N_r = 900 receptors located at the points

$$\begin{equation*}{\mathbf{x}}_{lq}=10\left(\mathrm{sin}\left(l\pi /31\right)\mathrm{cos}\left(2q\pi /30\right),\mathrm{sin}\left(l\pi /31\right)\mathrm{sin}\left(2q\pi /30\right),\mathrm{cos}\left(\pi /31\right)\right),\end{equation*}$$

with l, q = 1, ..., 30. As in examples 3 and 4, to avoid inverse crimes, we add 2% noise to the scattered data using formula (4.1).

We apply the RLA with the damped Gauss–Newton method at each frequency. We used the same stopping criteria for the Gauss–Newton method as in examples 3 and 4. We set the damping parameter α = 0.1 for the damped Gauss–Newton method at the initial frequency k = 0.5, and α = 0.1/(k||h||) for all other frequencies, where ||h|| is the 2-norm of the vector of coefficients of the update. Regarding the number of modes in the polynomial representing the update h, we set it to be N_p = ⌊2k⌋.

As we apply the RLA and increase the frequency, we note that oscillations are introduced in the reconstruction. These oscillations are an effect of the Gibbs phenomenon. They occur due to the limited number of modes used to recover the sharp edges of the obstacle. To mitigate the effect of the oscillations, we introduce an extra step in our reconstruction algorithm. After finding the domain update step h for the damped Gauss–Newton method, we apply a low-pass Gaussian filter as follows:

$$\begin{equation*}h\left(t\right)={h}_{0}^{\mathrm{c}}+\sum _{j=1}^{{N}_{p}}{\ell }_{j}\left({h}_{j}^{\mathrm{c}}\mathrm{cos}\left(2\pi j\left(t-0.5\right)\right)+{h}_{j}^{\mathrm{s}}\mathrm{sin}\left(2\pi j\left(t-0.5\right)\right)\right)\end{equation*}$$

where the filter constants are ${\ell }_{j}=\mathrm{exp}\left(-{\left(j/{N}_{p}\right)}^{2}/{\sigma }^{2}\right)$, j = 1, ..., N_p and σ² is a constant to define the damping of the filter. The application of this low-pass Gaussian filter follows a similar logic as in [6], which uses a curve smoother developed in [4].

In figures 10(a), (c) and (e), we present the reconstruction obtained at k = 20 using no filter, filter with σ² = 0.5, and filter with σ² = 0.1, respectively. In figures 10(b), (d) and (f) we present the cross section of the original obstacle and of the reconstructions at k = 6, k = 12 and k = 20 using no filter, filter with σ² = 0.5, and filter with σ² = 0.1, respectively.

Zoom In Zoom Out Reset image size

The results clearly show that although the general shape of the obstacle can be recovered without filtering, it is hard to recover the sharp features of the obstacle with high resolution. Using an appropriate filter produces a sharp reconstruction of the obstacle with very few oscillations, although determining the exact filtering parameter is not an obvious question and will be investigated in detail in a future work.