A direct imaging method for the half-space inverse scattering problem with phaseless data

We propose a direct imaging method based on the reverse time migration method for finding extended obstacles with phaseless total field data in the half space. We prove that the imaging resolution of the method is essentially the same as the imaging results using the scattering data with full phase information when the obstacle is far away from the surface of the half-space where the measurement is taken. Numerical experiments are included to illustrate the powerful imaging quality.

1. Introduction. In this paper, we study a direct imaging method for finding the support of an unknown extended obstacle embedded in the half space using the amplitude of the total field, which is measured on the boundary of half space far away from the obstacle. The algorithm does not require any a priori information of the physical properties of the obstacle such as penetrable or non-penetrable, and for non-penetrable obstacles, the type of boundary conditions on the boundary of the obstacle.
Let the penetrable obstacle occupy a bounded Lipschitz domain D ⊂ R 2 + = {(x 1 , x 2 ) T : x 1 ∈ R, x 2 > 0} with ν the unit outer normal to its boundary Γ D . We assume the incident wave is emitted by a point source located at x s on the surface Γ 0 = {(x 1 ; x 2 ) T : x 1 ∈ R, x 2 = 0}. Let N (x, y) be the Neumann Green function of the Hemholtz equation in the half space, (1) N (x, y) = Φ(x, y) + Φ(x, y ), where Φ(x, y) = i 4 H (1) 0 (k|x − y|) is the fundamental solution of the Hemholtz equation, and y = (y 1 , −y 2 ) T is the image point of y = (y 1 , y 2 ) T ∈ R 2 + . The total field is u(x, x s ) = N (x, x s ) + u s (x, x s ), where u s (x, x s ) is the solution of ∆u s (x, x s ) + k 2 n(x)u s (x, x s ) = −k 2 (n(x) − 1)N (x, x s ) in R 2 + , (2) ∂u s ∂x 2 (x, x s ) = 0 on Γ 0 , √ r ∂u s ∂r − iku s → 0 as r → ∞, r = |x|, where k > 0 is the wave number, n(x) ∈ L ∞ (R 2 + ) is a positive scalar function which is equal to 1 outside the scatter D. Condition (4) is the well known Sommerfeld radiation condition which guarantees the uniqueness of the solution. In this paper, by the scattering problem or scattering solution we always mean the solution satisfies the Sommerfeld radiation condition (4) .
In the diffractive optics imaging and radar imaging systems, measurement of the intensity of the total field is much easier and cheaper than the phase information of the field [13,14]. There are two types of approaches proposed in the literature to solve the inverse scattering problems with phaseless data. The first kind of approach is first to apply the phase retrieval algorithm to extract the phase information of the scattering field from the measurement of the intensity and then use the retrieved full field data in the classical imaging algorithms, see e.g. [15,2,24]. The second one is the iterative method which directly minimizes the difference between the received amplitude of the far field pattern and the synthesized phaseless data in the least square sense [19,17,18,1,5].
The reverse time migration (RTM) method, which consists of back-propagating the complex conjugated scattering field into the background medium and computing the crosscorrelation between the incident wave field and the backpropagated field to output the final imaging profile, is nowadays a standard imaging technique widely used in seismic imaging [3]. In [7,8,9], the RTM method for reconstructing extended targets using acoustic, electromagnetic and elastic waves at a fixed frequency is proposed and studied. Without using the small inclusion or geometrical optics assumption previously made in the literature, it is shown in [7,8,9] that the RTM imaging function is the total far field pattern of the scattering solution with the imaginary part of the fundamental solution as the incoming wave. This implies the RTM imaging function will always peak on the boundary of the scatterer.
We also remark that several direct sampling schemes for locating multiple multiscale acoustic scatterers are proposed in a very general and practical setting in [20,23], which are fast and robust against measurement noise. We also refer to [21,22] for direct sampling approaches to recover multiscale electromagnetic scatterers located in a homogeneous space.
In [11], the first attempt to apply direct imaging method using the idea of RTM method with phaseless data is made. It is shown that the proposed phaseless imaging function can achieve the same resolution as the RTM method applied to full phase data in [7]. The idea of developing direct imaging methods based on RTM has been extended to the electromagnetic waves in [12].
In this paper we study the phaseless imaging method in the half space since in practical applications, one may only be able to access the amplitude of the scattering data in the limited aperture. In [10], the RTM method for finding the extended obstacle embedded in the half space from scattering data with full phase information has been proposed and the resolution analysis is also proved under proper conditions. The half-space imaging function in [10] can be viewed as the frequency domain counterpart of the time domain imaging function previously proposed in the literature [27,26].
We denote by Ω the sampling domain in which the obstacle is sought.
where d > 0 is the aperture. Let u s (x, x s ) be the scattering solution of (2) − (3). The RTM imaging function proposed in [10] using full phase scattering data for imaging extended targets in the half space is given by, It is shown in [10] that this imaging function always peaks on the upper boundary of the obstacle.
In this paper, inspired by [11] we propose the following imaging function based on RTM for imaging obstacles with only intensity measurement |u(x r , x s )|: We will show in Theorem 4.1 below that Therefore the imaging resolution of our new phaseless RTM algorithm is essentially the same as the imaging results using the scattering data with the full phase information when the sources and measurements are placed far away from the scatterer. The rest of this paper is organized as follows. In section 2 we introduce our halfspace RTM algorithm for imaging the obstacle with phaseless data. In section 3 we introduce some preliminary results. In section 4 we consider the resolution of our algorithm for imaging penetrable obstacles. In section 5 we extend our theoretical results to non-penetrable obstacles with sound soft, sound hard, and impedance boundary condition. In section 6 we report several numerical examples to show the competitive performance of our phaseless RTM algorithm.

2.
Reverse time migration algorithm. In this section, we introduce the phaseless RTM method for inverse acoustic scattering problems in the half space. Assume that there are N s sources and N r receivers uniformly distributed on where d > 0 is the aperture. We always assume the obstacle D ⊂ Ω.
Let u i (x, x s ) = N (x, x s ) be the incident wave with source at x s ∈ Γ d 0 , and |u(x r , x s )| = |u s (x r , x s )+u i (x r , x s )| be the phaseless total data received at x r ∈ Γ d 0 , where u s (·, x s ) is the scattering solution of (2)-(3). We always assume that x s = x r for all s = 1,. . .,N s , r = 1,. . .,N r , to avoid the singularity of the incident field u i (x, x s ) at x = x r . This assumption can be easily satisfied in practical applications.
The phaseless RTM algorithm consists of back-propagating the corrected data ∆(x r , x s ) into the sampling domain using ∂Φ(xr,z) ∂x2(xr) , and then computing the imaginary part of the cross-correlation between ∂Φ(xs,z) ∂x2(xs) and the back-propagated field. Algorithm 2.1. (Phaseless RTM in the Half Space) Given the data |u(x r , x s )|=|u s (x r , x s ) + u i (x r , x s )| which is the measured amplitude of the total field at x r ∈ Γ d 0 , when the point source is emitted at x s ∈ Γ d 0 . 1 • Back-propagation: For s = 1, . . . , N s , compute the back-propagation field 2 • Cross-correlation: For z ∈ Ω, compute It is easy to see that This formula is used in all the numerical experiments in section 6. By letting N s , N r → ∞, we know that it can be viewed as an approximation of the following integral: We will show in section 4 that the above integral is indeed absolutely convergent.
3. Preliminary results. We first introduce some notations. For any bounded domain U ⊂ R 2 , let We introduce the following stability estimate of the forward acoustic scattering problem for penetrable obstacle in the half space, which can be proved by the same argument as in [4,Theorem 5.26].
has a unique radiating solution w ∈ H 1 (D). Moreover, there exists a constant C > 0, such that By Lemma 3.1, we can introduce the solution operator S : L 2 (D) → H 1 (D) by w = Sg which maps the incident wave g ∈ L 2 (D) to the scattering field w ∈ H 1 (D). In this paper we will denote by S the operator norm of S which depends generally on k, n(x) and the domain D.
For any y, z ∈ Ω, define Therefore, F (z, y) shares the same property as the point spread function: peaks when z = y and decays as the |z − y| becomes large. The following theorem for the imaging function I RTM (z) in (5) is given in [10, Theorem 5.2].
Theorem 3.2. For any z ∈ Ω, let ψ(y, z) be the radiation solution of the problem Then we have, for any z ∈ Ω, We remark that for the penetrable scatterers, ψ(y, z) is the scattering solution with the incoming field F (z, y). By (14) we expect the imaging function I RTM (z) will peak on the boundary of the scatterer and decay away from the scatterer. For the sound soft obstacles, when kh 1 and d h, it is shown in [10] by using the stationary phase theorem and Kirchhoff approximation that one cannot imaging the back part of the obstacle with the scattering data collected only on Γ 0 . The extensive numerical examples in [10] confirm the effectiveness of the RTM imaging function I RTM (z) in (5).

4.
Resolution analysis for the phaseless RTM algorithm. In this section, we prove the following theorem which shows that our RTM algorithm for the halfspace phaseless scattering data is asymptotically the same as the RTM algorithm using the scattering data with full phase information. In this section we make the following assumption: Here the constant C may dependent on kd D and n(·) L ∞ (D) , but is independent of k, h, d D .
We remark that S includes complicated dependence on k and the domain D in the stability estimate of the forward scattering problem. The constant C in the theorem may depend on the dimensionless constant kd D , but it does not depend on k or d D separately.
The proof of Theorem 4.1 depends on several lemmas that follow. We first observe that This yields Our goal now is to show that I 2 , I 3 are small when kh 1. The following estimates for the Hankel functions are proved in [6, (1.22 where d(x, D) = min y∈D |x − y| is the distance between x ∈ Γ d 0 and D, and the constant C may depend on kd D , n(·) L ∞ (D) , but is independent of k, d D .
Proof. By the integral representation theorem This implies by Lemma 3.1 that

By Lemma 4.2 we have for any
This completes the proof.
Thus by using the transform t = ( On the other hand, This completes the proof. Now we can give an estimate of the second term in the right hand side (16).
which implies by Lemma 4.4 that where we have used the fact that z 2 ≤ Ch for z ∈ Ω by the assumption (15).
Next we estimate the integral in we obtain by using Lemma 4.3 again that Thus by Lemma 4.4 we have This completes the proof by combining the estimates (18)- (19).
We now start at estimating the third term of the right hand side (16). Denote by δ = (h/k) 1 2 , and we introduce the set We will also use the notation that for x s , x r ∈ R 2 + , x s = (x 1s , x 2s ), x r = (x 1r , x 2r ).
Lemma 4.6. We have, for z ∈ Ω, where C may depend on kd D , n(·) L ∞ (D) , but is independent of k, h, d D .
Proof. It is obvious that Since |x s − z| ≥ z 2 for x s ∈ Γ d 0 , by (17), Lemma 4.3 and Lemma 4.4, we obtain where we have used δ = (h/k) 1/2 . The other terms can be estimated similarly. This completes the proof.
To proceed the estimate of the integral in Γ d 0 ×Γ d 0 \Q δ , we first recall the following von der Corput lemma for the oscillatory integral [16, P.152], [10, Lemma 3.2]. Lemma 4.7. Let −∞ < a < b < ∞, λ > 0, and u is a C 2 function in (a, b). If |u (t)| ≥ 1 for t ∈ (a, b) and u is monotone in (a, b), then for any φ defined in (a, b) where C may depend on kd D , n(·) L ∞ (D) , but is independent of k, h, d D , a, b.
where C may be dependent on kd D , n(·) L ∞ (D) , but is independent of k, h, d D , a, b.
Proof. We only prove the first estimate. The second one is similar. Let be the linear superposition of the scattering wave u s (·, x s ) along the source direction. Then w(y, z) is the scattering solution to (12) with This completes the proof.
where C may be dependent on kd D , n(·) L ∞ (D) , but is independent of k, h, d D .
Proof. Notice that u i (x r , x s ) = 2Φ(x r , x s ), by the asymptotic formula (21), we have

by Lemma 4.3 and Lemma 4.4 we have
We are now estimating the following integral By Lemma 4.9, we have

5.
Extensions. In this section we consider the reconstruction of sound-soft and impedance obstacles in the half space with phaseless data. For sound-soft obstacles, the measured data |u(x r , x s )| = |u s (x r , x s ) + N (x r , x s )|, where u s (x r , x s ) is the scattering solution of the following problem, By modifying the argument in Section 4 we can show the following theorem whose proof is omitted. Here S is the Dirichlet to Neumann map, and the constant C may depend on kd D , n(·) L ∞ (D) , but is independent of k, h, d D .
For impedance obstacles, the measured data |u(x r , x s )| = |u s (x r , x s )+N (x r , x s )|, where u s (x r , x s ) is the scattering solution of the following problem, By modifying the argument in Section 4 we can show the following theorem whose proof is omitted.
Here S is the Dirichlet to Neumann map, and the constant C may depend on kd D , n(·) L ∞ (D) , but is independent of k, h, d D .
For sound-soft or impedance obstacles, the studies in [10] show that the imaging function I RTM will peak at the boundary of the scatterer if kh 1 and d h. Therefore we again expect the imaging function I phaseless RTM (z) will have contrast on the boundary of the scatterer and decay outside the scatterer also for imaging impedance or sound-soft scatterers. 6. Numerical examples. In this section, we show several numerical examples by using the synthesized scattering data which are computed by standard Nyström methods. We use a uniform mesh with ten points per wavelength over the boundary to solve the boundary integral equation. The boundaries of the obstacles used in our numerical experiments are parameterized as follows: Circle : x 1 = ρ cos θ, x 2 = ρ sin θ, p − leaf : x 1 = r(θ) cos θ, x 2 = r(θ) sin θ, where r(θ) = 1 + 0.2 cos(pθ), P eanut : x 1 = cos θ + 0.2 cos 3θ, x 2 = sin θ + 0.2 sin 3θ, Rounded square : x 1 = cos 3 θ + cos θ, x 2 = sin 3 θ + sin θ.
In all our examples, we always choose h = 10 and d = 50. The sources x s and receivers x r are uniformly distributed on Γ d 0 with where s = 1, . . . , N s , r = 1, . . . , N r .
We consider the imaging of penetrable obstacles with different shapes including circle, peanut, kite and rounded square using phaseless data. Let the diffraction index n(x) = 1 outside D and n(x) = 0.5 inside D. The sampling domain is Ω = [−2, 2] × [8,12], with uniform mesh 201 × 201. The probe wave number is k = 4π, and N s = 512, N r = 512.
The top row of Figure 1 shows the imaging results of penetrable obstacles with different shapes. We observe that the location and upper boundary of the obstacle can be found effectively for different shapes. The interesting point is that the lower boundary of obstacles can also be imaged due to the diffraction energy back to the receivers after penetrating the media. As a comparison, we also compute the imaging results with the full phase scattering data which is shown in the bottom row of Figure 1. The imaging results with our method using phaseless data are almost the same as results using full phase scattering data, which confirms our theory in this paper.  Figure 2 shows that the location and upper boundary of the obstacles with different boundary conditions can be found without knowing any phase information of the data. Unlike the penetrable obstacle, the lower boundary of the obstacles can not be located as the scattering signals induced by the lower part of the obstacle propagate out without being recorded by receivers.
In our first test, we choose a single penetrable obstacle as the imaging object and the probe wavenumber is k = 4π for imaging with single frequency data. The top row of Figure 3 shows the imaging results with the noise level µ = 10%, 20%, 30%, 40% in the single frequency scattered phaseless data, respectively. The imaging results are still quite stable with the increasing of the noise level. We then generate multi-frequency data using the probe wavenumber k = 2π × [1 : 0.25 : 3]. As shown in the bottom row of Figure 3, the noise has been sucessfully suppressed and some undesired oscillations in imaging result with single frequency data are also cancelled by using multi-frequency data, which improves the imaging resolution of the imaging result.
For the second test, we choose two sound soft obstacles as imaging objects and the probe wavenumer is k = 4π for imaging with single frequency data. The search domain is Ω = [−5, 5]× [6,16] with a sampling mesh 201×201. The results shown in Figure 4 are imaged with phaseless data correlated with the nosie level as our previous test from left to right. Similarly, we also observe the noise can be suppressed and imaging results are improved by using multi-frequecy data.

Example 6.4.
We switch the position of two obstacles in Example 6.3 to show the impact of the limited data acquisition aperture. All the parameters are the same as in Example 6.3. The top row in Figure 5 are imaging results of two closely located obstacles with single frequency data (left plot) and multi-frequency data (right plot). We observe that the upper boundary of the larger obstacle is imaged as expected, but the smaller elliptic one below is completely invisible even we use the multi-frequency data. In the bottom row of Figure 5, we show the imaging results when these two obstacles are separated with a larger distance. We    can image the upper boundary of the elliptic obstacle below even only using single frequency data. The imaging result is also improved by using multi-frequency data. This indicates that the bottom obstacle can be imaged as long as the incident waves are not blocked by the larger top obstacle and the reflected wa ves can be acquisited on the surface.