Phaseless Imaging by Reverse Time Migration: Acoustic Waves This work is supported by National Basic Research Project under the grant 2011CB309700 and China NSF under the grants 11021101 and 11321061.

# Phaseless Imaging by Reverse Time Migration: Acoustic Waves ††thanks: This work is supported by National Basic Research Project under the grant 2011CB309700 and China NSF under the grants 11021101 and 11321061.

Zhiming Chen LSEC, Institute of Computational Mathematics and Scientific Engineering Computing, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, Beijing 100190, P.R. CHINA (zmchen@lsec.cc.ac.cn).    Guanghui Huang The Rice Inversion Project, Department of Computational and Applied Mathematics, Rice University, Houston, TX 77005-1892, USA (ghhuang@rice.edu). Previously: LSEC, Institute of Computational Mathematics and Scientific Engineering Computing, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, Beijing 100190, P.R. CHINA (ghhuang@lsec.cc.ac.cn).
###### Abstract

We propose a reliable direct imaging method based on the reverse time migration for finding extended obstacles with phaseless total field data. 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. The imaginary part of the cross-correlation imaging functional always peaks on the boundary of the obstacle. Numerical experiments are included to illustrate the powerful imaging quality.

## 1 Introduction

We consider in this paper inverse scattering problems with phaseless data which aim to find the support of unknown obstacles embedded in a known background medium from the knowledge of the amplitude of the total field measured on a given surface far away from the obstacles. Let the sound soft obstacle occupy a bounded Lipschitz domain with the unit outer normal to its boundary . Let be the incident wave and the total field is with being the solution of the following acoustic scattering problem:

 Δus+k2us=0    in R2∖¯D, \hb@xt@.01(1.1) us=−ui     on ΓD, \hb@xt@.01(1.2) √r(∂us∂r−ikus)→0  as   r=|x|→+∞, \hb@xt@.01(1.3)

where is the wave number. The condition (LABEL:p3) is the outgoing Sommerfeld radiation condition which guarentees the uniqueness of the solution. In this paper, by the radiation or scattering solution we always mean the solution satisfies the Sommerfeld radiation condition (LABEL:p3). For the sake of the simplicity, we mainly consider the imaging of sound soft obstacles in this paper. Our algorithm does not require any a priori information of the physical properties of the obstacles such as penetrable or non-penetrable, and for non-penetrable obstacles, the type of boundary conditions on the boundary of the obstacle. The extension of our theoretical results for imaging other types of obstacles will be briefly considered in section 4.

In the diffractive optics imaging and radar imaging systems, it is much easier to measure the intensity of the total field than the phase information of the field [9, 10, 26]. It is thus very desirable to develop reliable numerical methods for reconstructing obstacles with only phaseless data, that is, the amplitude information of the total field . In recent years, there have been considerable efforts in the literature to solve the inverse scattering problems with phaseless data. One approach is to image the object with the phaseless data directly in the inversion algorithm, see e.g. [10, 19]. The other 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. [11]. We also refer to [1] for the continuation method and [17, 13, 14] for inverse scattering problems with the data of the amplitude of the far field pattern. In [15] some uniqueness results for phaseless inverse scattering problems have been obtained.

The reverse time migration (RTM) method, which consists of back-propagating the complex conjugated scattering field into the background medium and computing the cross-correlation 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 [2]. In [3, 4, 5], the RTM method for reconstructing extended targets using acoustic, electromagnetic and elastic waves at a fixed frequency is proposed and studied. The resolution analysis in [3, 4, 5] is achieved without using the small inclusion or geometrical optics assumption previously made in the literature.

In this paper we propose a direct imaging algorithm based on reverse time migration for imaging obstacles with only intensity measurement with point source excitations. We prove that the imaging resolution of the new algorithm is essentially the same as the imaging results using the scattering data with the full phase information, that is, our imaging function always peaks on the boundary of the obstacles. To the best knowledge of the authors, our method seems to be the first attempt in applying non-iterative method for reconstructing obstacles with phaseless data except [9] in which a direct method is considered for imaging a penetrable obstacle under Born approximation using plane wave incidences. We will extend the RTM method studied in this paper for electromagnetic probe waves in a future paper.

The rest of this paper is organized as follows. In section 2 we introduce our RTM algorithm for imaging the obstacle with phaseless data. In section 3 we consider the resolution of our algorithm for imaging sound soft obstacles. In section 4 we extend our theoretical results to non-penetrable obstacles with the impedance boundary condition and penetrable obstacles. In section 5 we report several numerical experiments to show the competitive performance of our phaseless RTM algorithm.

## 2 Reverse time migration method

In this section we introduce the RTM method for inverse scattering problems with phaseless data. Assume that there are emitters and receivers uniformly distributed respectively on and , where are the disks of radius respectively. We denote by the sampling domain in which the obstacle is sought. We assume the obstacle and is inside .

Let , where is the fundamental solution of the Helmholtz equation with the source at , be the incident wave and be the phaseless data received at , where is the solution to the problem (LABEL:p1)-(LABEL:p3) with . We additionally assume that for all , to avoid the singularity of the incident field at . This assumption can be easily satisfied in practical applications. In the following, without loss of generality, we assume .

Our RTM algorithm consists of back-propagating the corrected data:

 Δ(xr,xs)=|u(xr,xs)|2−|ui(xr,xs)|2ui(xr,xs) \hb@xt@.01(2.1)

into the domain using the fundamental solution and then computing the imaginary part of the cross-correlation between and the back-propagated field.

###### Algorithm 2.1

(RTM for Phaseless data) Given the data which is the measurement of the total field at when the point source is emitted at , , .
Back-propagation: For , compute the back-propagation field

 vb(z,xs)=−2πRrNrNr∑r=1Φ(xr,z)Δ(xr,xs),    ∀ z∈Ω. \hb@xt@.01(2.2)

Cross-correlation: For , compute

 I(z)=−k2Im{2πRsNsNs∑s=1ui(z,xs)vb(z,xs)}. \hb@xt@.01(2.3)

It is easy to see that

 I(z)=−k2Im{(2π)2RsRrNsNrNs∑s=1Nr∑r=1Φ(z,xs)Φ(xr,z)Δ(xr,xs)},  ∀ z∈Ω. \hb@xt@.01(2.4)

This is the formula used in our numerical experiments in section 5. By letting , we know that (LABEL:scor2) can be viewed as an approximation of the following continuous integral:

 ^I(z)=−k2Im∫Γs∫ΓrΦ(z,xs)Φ(xr,z)Δ(xr,xs)ds(xr)ds(xs),  ∀ z∈Ω. \hb@xt@.01(2.5)

We remark that the above RTM imaging algorithm is the same as the RTM method in [3] except that the input data now is instead of . Hence, the code of the RTM algorithm for imaging the obstacle with phaseless data requires only one line change from the code of the RTM method for imaging the obstacle with full phase information.

## 3 The resolution analysis

In this section we study the resolution of the Algorithm 2.1. We first introduce some notation. For any bounded domain with Lipschitz boundary , let be the weighted norm and be the weighted norm, where is the diameter of and

 |v|12,Γ=(∫Γ∫Γ|v(x)−v(y)|2|x−y|2ds(x)ds(y))1/2.

By scaling argument and trace theorem we know that there exists a constant independent of such that for any ,

 ∥ϕ∥H1/2(ΓD)+∥∂ϕ/∂ν∥H−1/2(ΓD)≤Cmaxx∈¯D(|ϕ(x)|+dD|∇ϕ(x)|). \hb@xt@.01(3.1)

The following stability estimate of the forward acoustic scattering problem is well-known [8, 20].

###### Lemma 3.1

Let , then the scattering problem:

 Δw+k2w=0in R2∖¯D,    w=g    on ΓD, \hb@xt@.01(3.2) √r(∂w∂r−ikw)→0,    as r→∞, \hb@xt@.01(3.3)

admits a unique solution . Moreover, there exists a constant such that .

The far field pattern , where , of the solution of the scattering problem (LABEL:ha)-(LABEL:ha1) is defined as (cf. e.g. [8, P. 67]):

 w∞(^x)=eiπ4√8πk∫ΓD[w(y)∂e−ik^x⋅y∂ν(y)−∂w(y)∂ν(y)e−ik^x⋅y]ds(y). \hb@xt@.01(3.4)

It is well-known that for the scattering solution of (LABEL:ha)-(LABEL:ha1) (cf. e.g. [3, Lemma 3.3])

 −Im∫ΓDw∂¯w∂νds=k∫S1|w∞(^x)|2d^x. \hb@xt@.01(3.5)

Now we turn to the analysis of the imaging function in (LABEL:scor3). We first observe that

 Δ(xr,xs)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯us(xr,xs)+|us(xr,xs)|2ui(xr,xs)+us(xr,xs)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ui(xr,xs)ui(xr,xs). \hb@xt@.01(3.6)

This yields

 ^I(z) = −k2Im∫Γs∫ΓrΦ(z,xs)Φ(xr,z)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯us(xr,xs)ds(xr)ds(xs) \hb@xt@.01(3.7) −k2Im∫Γs∫ΓrΦ(z,xs)Φ(xr,z)|us(xr,xs)|2ui(xr,xs)ds(xr)ds(xs) −k2Im∫Γs∫ΓrΦ(z,xs)Φ(xr,z)us(xr,xs)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ui(xr,xs)ui(xr,xs)ds(xr)ds(xs).

The first term is the RTM imaging function with full phase information in [3] and thus can be analyzed by the argument there. Our goal now is to show the last two terms at the right hand side of (LABEL:c3) are small. We start with the following lemma.

###### Lemma 3.2

We have for any .

Proof. We first recall the following estimates for Hankel functions [6, (1.22)-(1.23)]:

 |H(1)0(t)|≤(2πt)1/2,  |H(1)1(t)|≤(2πt)1/2+2πt,    ∀t>0. \hb@xt@.01(3.8)

By the integral representation formula, we have

 us(xr,xs)=∫ΓD(us(y,xs)∂Φ(xr,y)∂ν(y)−∂us(y,xs)∂ν(y)Φ(xr,y))ds(y). \hb@xt@.01(3.9)

By (LABEL:d00) we have

 ∥Φ(xr,⋅)∥H1/2(ΓD)+∥∂Φ(xr,⋅)/∂ν∥H−1/2(ΓD)≤C(1+kdD)(kRr)−1/2.

The lemma follows now from Lemma LABEL:lem:wp and the fact that for .

###### Lemma 3.3

We have for any .

Proof. We use the following integral formula [22, 6]

 H(1)0(t)=−2iπeit∫∞0e−rtr1/2(r−2i)1/2dr,    ∀t>0,

where for . By the change of variable

 |H(1)0(t)|≥2πRe∫∞0e−ss1/2(s−2it)1/2ds = 2π∫∞0e−ss1/2|s−2it|√|s−2it|+s2ds ≥ 2π∫∞0e−s|s−2it|ds ≥ 25π∫1te−ssds,

where in the last inequality we have used for . This completes the proof by noticing that .

The following lemma gives an estimate of the second term at the right hand side of (LABEL:c3).

###### Lemma 3.4

We have

 ∣∣ ∣∣k2∫Γs∫ΓrΦ(z,xs)Φ(xr,z)|us(xr,xs)|2ui(xr,xs)ds(xr)ds(xs)∣∣ ∣∣≤C(1+kdD)4(kRs)−1/2.

Proof. Denote . Let , , where . Since , , we obtain easily

 |xr−xs|=Rs√1+τ2−2τcos|θr−θs|≥2Rs√τsin|θr−θs|2 \hb@xt@.01(3.10)

by Cauchy-Schwarz inequality. Now for , we have then either or , where .

By (LABEL:b1) and Lemma LABEL:lem:4.1,

 ∣∣ ∣∣∫∫ΩkΦ(z,xs)Φ(xr,z)|us(xr,xs)|2ui(xr,xs)ds(xr)ds(xs)∣∣ ∣∣ \hb@xt@.01(3.11) ≤ C(1+kdD)4(kRs)−3/2(kRr)−3/2∫∫Ωk|ln(k|xr−xs|)|ds(xr)ds(xs).

By (LABEL:y2) we have

 ∫∫Ωk|ln(k|xr−xs|)|ds(xr)ds(xs) ≤ −∫∫Ωkln(2kRs√τsin|θr−θs|2)dθrdθs ≤ 2πRrRs∣∣∣∫θ00ln(2kRs√τsint)dt∣∣∣ ≤ CRrRs,

where we have used integration by parts in obtaining the last inequality. Substituting the above estimate into (LABEL:c4) we obtain

 ∣∣ ∣∣∫∫ΩkΦ(z,xs)Φ(xr,z)|us(xr,xs)|2ui(xr,xs)ds(xr)ds(xs)∣∣ ∣∣≤Ck−2(1+kdD)4(kRs)−1. \hb@xt@.01(3.12)

Next we estimate the integral in . Since is an increasing function of [24, p. 446], we have for , , and thus

 |xr−xs||ui(xr,xs)|2≥132k−1∣∣∣H(1)0(12)∣∣∣2=Ck−1,

which implies by using Lemma LABEL:lem:4.1 and (LABEL:b1) again that

 ∣∣ ∣∣∫∫Γr×Γs∖¯ΩkΦ(z,xs)Φ(xr,z)|us(xr,xs)|2ui(xr,xs)ds(xr)ds(xs)∣∣ ∣∣≤Ck−2(1+kdD)4(kRs)−1/2. \hb@xt@.01(3.13)

This completes the proof by combining the above estimate with (LABEL:c5).

Now we turn to the estimation of the third term at the right hand side of (LABEL:c3). Denote and and .

###### Lemma 3.5

We have

 ∣∣ ∣∣k2∫∫QδΦ(z,xs)Φ(xr,z)us(xr,xs)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ui(xr,xs)ui(xr,xs)ds(xr)ds(xs)∣∣ ∣∣≤C(1+kdD)2(kRs)−1/2.

Proof. The proof follows easily from Lemma LABEL:lem:4.1 and (LABEL:b1) and the fact that .

To estimate the integral in , we recall first the following useful mixed reciprocity relation [16], [23, P.40].

###### Lemma 3.6

Let . Then for any , where is the far field pattern in the direction of the scattering solution of (LABEL:p1)-(LABEL:p3) with and is the scattering solution of (LABEL:p1)-(LABEL:p3) with the incident plane wave .

###### Lemma 3.7

Let be the far field pattern in the direction of the scattering solution of the Helmholtz equation with the incident plane wave . Then .

Proof. The proof follows from the definition of the far field pattern in (LABEL:far) with on , Lemma LABEL:lem:wp, and (LABEL:d00). Here we omit the details.

The following lemma is essentially proved in [8, Theorem 2.5].

###### Lemma 3.8

For any , the solution of the scattering problem (LABEL:ha)-(LABEL:ha1) satisfies the asymptotic behavior:

 w(x)=eik|x|√|x|w∞(^x)+γ(x),

where .

Proof. First we have the following integral representation formula

 w(x)=∫ΓD[w(y)∂Φ(x,y)∂ν(y)−Φ(x,y)∂w(y)∂ν(y)]ds(y),    ∀x∈R2∖¯D. \hb@xt@.01(3.14)

By the asymptotic formulae of Hankel functions [24, P.197], for ,

 H(1)n(t)=(2πt)1/2ei(kt−nπ2−π4)+Rn(t),  |Rn(t)|≤Ct−3/2,  ∀t>0, \hb@xt@.01(3.15)

and the simple estimate for any , , we have

 Φ(x,y) = \hb@xt@.01(3.16) ∂Φ(x,y)∂ν(y) = \hb@xt@.01(3.17)

where for some constant independent of and . The proof completes by inserting (LABEL:r3)-(LABEL:r4) into (LABEL:r1) and using Lemma LABEL:lem:wp and (LABEL:d00). Here we omit the details.

We also need the following slight generalization of Van der Corput lemma for the oscillatory integral [12, P.152].

###### Lemma 3.9

For any , for every real-valued function that satisfies for . Assume that is a division of such that is monotone in each interval , . Then for any function defined on with integrable derivative, and for any ,

 ∣∣∣∫baeiλu(t)ϕ(t)dt∣∣∣≤(2N+2)λ−1[|ϕ(b)|+∫ba|ϕ′(t)|dt].

Proof. By integration by parts we have

 ∫baeiλu(t)dt=[eiλu(t)iλu′(t)]ba−∫baeiλu(t)ddt(1iλu′(t))dt.

Since is monotone in each interval , , and in , we have

 ∣∣ ∣∣∫baeiλu(t)ddt(1iλu′(t))dt∣∣ ∣∣≤N∑i=1λ−1∣∣ ∣∣∫xixi−1ddt(1u′(t))dt∣∣ ∣∣≤2Nλ−1,

which implies . For the general case, we denote and use integration by parts to obtain

 ∫baϕ(t)eiλu(t)dt=ϕ(b)F(b)−∫baF(t)ϕ′(t)dt.

This completes the proof by using .

###### Lemma 3.10

We have

 ∣∣ ∣∣k2∫∫Γr×Γs∖¯QδΦ(z,xs)Φ(xr,z)us(xr,xs)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ui(xr,xs)ui(xr,xs)ds(xr)ds(xs)∣∣ ∣∣ ≤C(1+kdD)3(kRs)−1/2+C(1+k|z|)2(kRs)−1.

Proof. We first observe that for , we have

 k|xr−xs|≥2kRs√τ∣∣∣sinθr−θs2∣∣∣≥2kRs√τsinδ2≥12(kRs)1/2√τ,

where we have used the fact that for . Thus by (LABEL:r2) we obtain

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ui(xr,xs)ui(xr,xs)=e−2ik|xr−xs|+iπ2+ρ0(xr,xs), \hb@xt@.01(3.18)

where . Similar to (LABEL:r3) we have

 Φ(z,xs) = \hb@xt@.01(3.19) Φ(z,xr) = \hb@xt@.01(3.20)

where , . Next by Lemma LABEL:lem:far, the mixed reciprocity in Lemma LABEL:lem:4.5, and (LABEL:d00) we have

 us(xr,xs) = eikRr√Rru∞ps(^xr,xs)+ρ2(xr,xs) \hb@xt@.01(3.21) = eikRr√Rrγmus(xs,−^xr)+ρ2(xr,xs) = eik(Rr+Rs)√RrRsγmu∞pl(^xs,−^xr)+ρ2(xr,xs)+ρ3(xr,xs),

where is the far field pattern of the scattering solution of the Helmholtz equation with the incident plane wave and

 |ρ2(xr,xs)| ≤ C(1+kdD)3(kRr)−3/2∥Φ(⋅,xs)∥H1/2(ΓD) ≤ C(1+kdD)4(kRr)−3/2(kRs)−1/2, |ρ3(xr,xs)| ≤ C(1+kdD)3(kRr)−1/2(kRs)−3/2∥e