A reference ball based iterative algorithm for imaging acoustic obstacle from phaseless far-field data

In this paper, we consider the inverse problem of determining the location and the shape of a sound-soft obstacle from the modulus of the far-field data for a single incident plane wave. By adding a reference ball artificially to the inverse scattering system, we propose a system of nonlinear integral equations based iterative scheme to reconstruct both the location and the shape of the obstacle. The reference ball technique causes few extra computational costs, but breaks the translation invariance and brings information about the location of the obstacle. Several validating numerical examples are provided to illustrate the effectiveness and robustness of the proposed inversion algorithm.


Introduction
It is well known that inverse scattering problems (ISP) play a significantly crucial role in a vast variety of realistic applications such as radar sensing, ultrasound tomography, biomedical imaging, geophysical exploration and noninvasive detecting. Therefore, effective and efficient numerical inversion approaches have been extensively and intensively studied in the recent decades [4].
In typical inverse scattering problems, reconstruction of the geometrical information is based on the knowledge of full measured data (both the intensity and phase information are collected). However, in a number of practical scenarios, to measure the full data is extremely difficult or might even be unavailable. Thus, phaseless inverse scattering problems arise naturally and attract a great attention from mathematical and numerical point of view. Without the phase information, the theoretical justifications of the uniqueness issue and the development of inversion algorithms are more challenging than the phased inverse scattering problem.
Various numerical methods have been proposed for solving phaseless inverse acoustic obstacle scattering problems. Kress and Rundell [13] investigate the phaseless inverse obstacle scattering and propose a Newton method for imaging a two-dimensional sound-soft obstacle from the modulus of the far-field data with only one incoming direction. In particular, it is pointed out that the location of the obstacle cannot be reconstructed since the modulus of the far-field pattern has translation invariance. This means that the solution of the inverse problem is not unique. If the location of the obstacle is known, then a uniqueness result is presented in [20] for a special case, that is, the sound-soft ball can be determined by phaseless far-field pattern when the product of wavenumber and the radius of the ball is less than a constant. In [8], a nonlinear integral equation method is investigated for the two-dimensional shape reconstruction from a given incident field and the modulus of the far-field data. This method involves the full linearization of a pair of integral equations, i.e., field equation and phaseless data equation, with respect to both the boundary parameterization and the density. Then, the nonlinear integral equation method is extended to the three-dimensional shape reconstruction in [10]. The problem is divided into two subproblems including the shape reconstruction from the modulus of the farfield data and the location identification based on the invariance relation using only few far field measurements in the backscattering direction. In addition, fundamental solution method [11] and a hybrid method [17] are proposed to detect the shape of a sound-soft obstacle by using of the modulus of the farfield data for one incident field.
Besides the aforementioned phaseless inverse acoustic obstacle scattering, there exist other types of phaseless inverse scattering problems as well as the relevant numerical methods. Based on the physical optics approximation and the local maximum behavior of the backscattering far-field pattern, a numerical method is developed in [19] to reconstruct an electromagnetic polyhedral PEC obstacle from a few phaseless backscattering measurements. In [5], a direct imaging method based on reverse time migration is proposed for reconstructing extended obstacles by phaseless electromagnetic scattering data. A continuation method [3] is developed to reconstruct the periodic grating profile from phaseless scattered data. In [2], a recursive linearization algorithm is proposed to image the multi-scale rough surface with tapered incident waves from measurements of the multi-frequency phaseless scattered data. In [24], an iterative method based on Rayleigh expansion is developed to solve the inverse diffraction grating problem with super-resolution by using phase or phaseless near-field data. For a recent work on the Fourier method for solving phaseless inverse source problem, we refer to [23]. We also refer to [1,14,15,16] for phaseless inverse medium scattering problems.
A major difficulty in solving the phaseless inverse obstacle scattering problem is to determine the location of the obstacle from phaseless farfield data, since the modulus of the far-field pattern is invariant under translations. Recently, a recursive Newton-type iteration method is developed in [21] to recover both the location and the shape of the obstacle simultaneously from multi-frequency phaseless far-field data. In this numerical method, some sets of superpositions of two plane waves with different directions are used as the incident fields to break the translation invariance property of the phaseless far-field pattern.
In this paper, we consider the inverse problem of determining the location and the shape of a sound-soft obstacle from the modulus of the far-field data for a single incident plane wave. In a recent work [6], the nonlinear integral equation method proposed by Johansson and Sleeman [7] is extended to reconstruct the shape of a sound-soft crack by using phaseless far-field data for one incident plane wave. Motivated by [6,7] and the reference ball technique in [18,22], we first introduce an artificial reference ball to the inverse scattering system, and propose an iterative scheme which involves a system of nonlinear and ill-posed integral equations to reconstruct both the location and the shape of the obstacle. Since the location of the reference ball is known and fixed, it has the capability of calibrating the scattering system so that the translation invariance does not occur. As a result, the location information of the obstacle could be recovered with negligible additional computational costs.
To our best knowledge, this is the first attempt in the literature towards reconstructing both the location and the shape of the unknown obstacle by using the intensity-only far-field data due to a single incident plane wave. The novelty of this work lies in the incorporation of the reference ball technique into the iterative scheme. Hence, our approach exhibits the following salient features: First, rather than the superposition of different incident waves, only a single incident plane wave with a fixed wavenumber and a fixed incoming direction is needed. Second, the iteration procedure is very fast and easy to implement because the Fréchet derivatives can be explicitly formulated and thus do not require the forward solver. Third, the location and profile information of the obstacle could be simultaneously reconstructed. Finally, the iterative scheme is robust in the sense that it is insensitive to the parameters such as the initial guess, the location and size of the reference ball, as well as the measurement noise.
The rest of this paper is organized as follows. In the next section, we formulate the phaseless inverse obstacle scattering problem in conjunction with the reference ball. In Section 3, we derive a system of boundary integral equations with a reference ball, and an iterative scheme is presented to solve the boundary integral equations. In Section 4, several numerical examples will be presented, including the numerical implementation details. Finally, some concluding remarks are summarized in Section 5.

Problem formulation
Let D ⊂ R 2 be an open bounded domain with C 2 boundary and the positive constant κ be the wavenumber. Given an incident plane wave u i = e iκx·d with the incoming direction d, the forward/direct obstacle scattering problem is to find the total field u satisfying the Helmholtz equation and the Dirichlet boundary condition The total field u = u i + u s is given as the sum of the known incident wave u i and the unknown scattered wave u s which is required to fulfill the Sommerfeld radiation condition Further, the scattered field u s has the following asymptotic behavior [4] uniformly for all directionx := x/|x|, where the complex-valued function u ∞ (x) defined on the unit circle Ω is known as the far field pattern or scattering amplitude. Then, the phaseless inverse scattering problem is stated as follows: Problem 1 (Phaseless ISP) Given an incident plane wave u i for a single wavenumber and a single incident direction d, together with the corresponding phaseless far-field data |u ∞ D (x)|,x ∈ Ω, due to the unknown obstacle D, determine the location and shape of ∂D.
It has been pointed out in [13] that Problem 1 does not admit a unique solution, namely, the location of the obstacle cannot be reconstructed since the modulus of the far-field pattern has translation invariance. Specifically, for the shifted domain D h := {x + h : x ∈ D} with a fixed vector h ∈ R 2 , the far-field pattern u ∞ D h satisfies the relation Moreover, the ambiguity induced by this translation invariance relation cannot be remedied by using a finite number of incident waves with different wavenumbers or different incident directions, see [13]. Our goal in this paper is to overcome this difficulty via the reference ball based iterative scheme. To this end, we reformulate Problem 1 as follows: Problem 2 (Phaseless ISP with a reference ball) Let B ⊂ R 2 be the artificially added sound-soft ball such that D ∩ B = ∅. Given an incident plane wave u i for a single wavenumber and a single incident direction d, together with the corresponding phaseless far-field data |u ∞ D∪B (x)|,x ∈ Ω, due to the scatterer D ∪ B, determine the location and shape of ∂D.
The geometry setting of Problem 2 is illustrated in Fig. 1. For brevity, we denote by Γ 1 := ∂D and Γ 2 := ∂B in what follows. In the next section, we will introduce a system of nonlinear integral equations based iterative scheme for solving Problem 2. After that, several numerical experiments will be conducted to demonstrate the feasibility of breaking the translation invariance using the proposed method.

Nonlinear integral equations
We begin this section with establishing a system of nonlinear integral equations for the inverse scattering problem. Denote the fundamental solution to the Helmholtz equation by where H 0 denotes the zero-order Hankel function of the first kind. By applying the Huygens' principle [4, Theorem 3.14] to the scattering system D ∪ B, we obtain Then the far-field pattern of the scattered field u s is given by where γ = e iπ/4 / √ 8κπ and ν denotes the unit outward normal. We define the single-layer operator (S jl g)(x) := Γj Φ(x, y)g(y) ds(y), x ∈ Γ l , j, l = 1, 2 and the corresponding far-field operator Let x ∈ R 2 \D ∪ B tend to the boundary Γ 1 and Γ 2 , respectively, and in terms of (4), the Dirichlet boundary condition (2) and the continuity of the single-layer potential, one can readily deduce the following field equations with respect to the densities g j = ∂u/∂ν| Γj , j = 1, 2.
On the other hand, let |x| → ∞, then the linearity of the forward scattering problem and equation (5) lead to the phaseless data equation

The iterative procedure
We now seek a sequence of approximation to Γ 1 by solving the field equations (6)- (7) and phaseless data equation (8) in an alternating manner. Given an approximation for the boundary Γ 1 one can solve (6)- (7) for g 1 and g 2 . Then keeping g 1 and g 2 fixed, the equation (8) is linearized with respect to Γ 1 to update the boundary approximation.
For the sake of simplicity, the boundary Γ 1 is assumed to be a starlike curve with the parametrized form the boundary Γ 2 and the unit circle Ω are parameterized by Let p j (t) be the points on Γ j , described by p 1 (t) := (c 1 , c 2 ) + r(t)(cos t, sin t) and p 2 (t) : . Further, we introduce the parameterized single-layer operators S jl and far field operator S ∞ j by and where we have set ψ j (τ ) = G j (τ )g j (p j (τ )). Here, G 1 (τ ) = r 2 (τ )+ dr dτ (τ ) 2 1/2 and G 2 (τ ) = R denote the Jacobian of the transformation. The corresponding right-hand sides are w l (t) = u i (p l (t)) and w ∞ (t) = u ∞ D∪B (x(t)). Thus we can obtain the parametrized integral equations (6)-(8) in the form The linearization of equation (14) with respect to p 1 requires the Fréchet derivatives of the operator A ∞ 1 , that is where the update q(τ ) = (∆c 1 , ∆c 2 ) + ∆r(τ )(cos τ, sin τ ). Hence, by using of the product rule, we have and the linearization of (14) leads to where As usual for iterative algorithms, the stopping criteria is necessary to justify the convergence in numerics. Regarding our iterative procedure, we choose the following relative error estimator for some sufficiently small parameter ǫ > 0 depending on the noise level. Here, p 1 is the kth approximation of the boundary Γ 1 .
We are now in the position to present the reference ball based iteration algorithm: Algorithm: Iterative procedure for phaseless inverse scattering Step 1 With a mild a priori information of the unknown scatterer D, add a suitable reference ball B such that D ∩ B = ∅; Step 2 Emanate an incident plane wave with a fixed wavenumber κ > 0 and a fixed incident direction d ∈ Ω, and then collect the corresponding noisy phaseless far-field data |u ∞ D∪B (x)|,x ∈ Ω for the scatterer D ∪ B; Step 3 Select an initial star-like curve Γ (0) for the boundary ∂D and the error tolerance ǫ. Set k = 0; Step 4 For the curve Γ (k) , find the densities ψ 1 and ψ 2 from (12)-(13).

Remark 1
The advantage of introducing the reference ball is that we would be able to obtain the location of the obstacle via this iterative method, since the update information (∆c 1 , ∆c 2 ) about the location of the obstacle is contained in the term ℜ(A ∞ 2 (p 2 , ψ 2 )A ′∞ 1 [p 1 , ψ 1 ]q) of equation (16).
Remark 2 For typical nonlinear integral equations based iterative schemes (see, e.g., [8]), the linearization are carried out with respect to both the boundary and the density functions, so the process of solving unknowns should be essentially intertwined. In contrast, we would like to emphasize that our inversion scheme requires the linearization only with respect to the boundary for phaseless data equation, hence the field equations and the phaseless data equation could be completely decoupled and then solved alternatively and separately in each iterative loop. Moreover, the Fréchet derivatives in this paper also does not rely on any solution process and can be explicitly given. Therefore, the numerical implementation is now significantly simplified.
In the numerical examples, we obtain the update ξ from a scaled Newton step with Tikhonov regularization and H 2 penalty term, that is where the scaling factor ρ ≥ 0 is fixed throughout the iterations. According to [9], the regularization parameters λ in equation (22) are chosen as Note that in each step of the iteration, the derivative dr/dτ is calculated by resorting to the approximation (18) and r (k) (τ ) = r (k−1) (τ ) + ∆r (k−1) (τ ), k = 1, 2 · · · . Analogously to [13], the initial approximation is chosen as a circle and the values of the Fourier modes in the directions cos θ and sin θ (i.e., coefficients α 1 and β 1 ) are set to be the exact values in the reconstructions. The reason for doing so is to ease the comparison with the exact curve. So the update procedure does not take these two modes into account.
In the subsequent figures, the exact boundary curves are displayed as solid lines, the reconstruction are depicted with the dashed lines −−, the initial guess are taken to be circles with radius r (0) = 0.1 indicated by the dash-dotted lines ·−. The incident directions are denoted by arrows. Throughout all the numerical examples, we set the wavenumber κ = 2, the scaling factor ρ = 0.6, and the parameter M = 5.
Here the incident direction d = (cos(−π/6), sin(−π/6)) is used in Fig. 2-4. Several snapshots of the iterative process are shown in Fig. 2 and Fig. 3, where the initial guess is the circle centered at (−0.7, 0.45) with radius 0.1, and     the reference ball is centered at (4, 0) with radius 0.4. Moreover, the relative L 2 error Er k between the reconstructed and exact boundaries and the error E k defined in (17) are also presented with respect to the number of iterations.
As we can see from the figures, the trend of two error curves is basically the same for larger number of iteration. So, the choice of the stopping criteria is reasonable. The reconstructions with different reference balls for the incoming wave direction d = (cos(−π/6), sin(−π/6)) and d = (cos(4π/3), sin(4π/3)) are shown in Fig. 4 and Fig. 6 respectively, and the reconstructions with different incoming wave directions are presented in Fig. 5. As shown in these results, the location and shape of the obstacle could be simultaneously and satisfactorily reconstructed.
In this example, the reconstructions with 1% noise and 5% noise, corresponding to the incoming wave direction d = (cos(2π/3), sin(2π/3)), are shown in Fig. 7 and Fig. 8, respectively. Here we use the initial guess (c 2 ) = (0.3, −0.6), r (0) = 0.1 and the reference ball (b 1 , b 2 ) = (4, 0), R = 0.4. The relative L 2 error Er k and the error E k are also presented in the figures. The reconstructions with different reference balls for the incoming wave direction d = (cos(π/6), sin(π/6)) are shown in Fig. 9, and the reconstructions with different initial guesses for the incoming wave direction d = (cos(2π/3), sin(2π/3)) are shown in Fig. 10. In this example, the reconstructed obstacle and the relative error Er k and E k from the phaseless far-field data with 1% noise and 5% noise, corresponding to the incoming wave direction d = (cos(π/6), sin(π/6)), are shown in Fig. 11. The relative L 2 error Er k and the error E k are also presented in the figures. The influences of the choices of initial guesses and reference balls are shown in Fig. 12 and Fig. 13 respectively.
The above numerical results illustrate that by adding a reference ball artificially to the inverse scattering system the iteration method gives a feasible reconstruction of the location and the shape of the obstacle from phaseless farfield data for one incident field. The reference ball causes few extra computational costs, but breaks the translation invariance and brings information about the location of the obstacle. In addition, a promising feature of the algorithm is that the Fréchet derivatives involved in our method can be explicitly characterized as integral operators and thus easily evaluated. Hence it does not require the solution of the forward problem in each iteration step and it is very easy to implement with computational efficiency. To evaluate the computational time, our codes of the experiments are written in Matlab and run on a laptop with 2.6 GHz CPU. According to the fact that the CPU time for all the reconstructions is less then 10 seconds, we conclude that our algorithm is very fast.

Conclusions and future works
In this paper, a new numerical method is devised to solve the inverse obstacle scattering problem from the modulus of the far-field data for one incident field. That is, by introducing a reference ball artificially to the inverse scattering system, the translation invariance of the phaseless far-field pattern can be broken down, and an iterative scheme which is based on a system of boundary integral equations is propose to reconstruct both the location and shape of the obstacle. The reference ball causes few extra computational costs, but breaks the translation invariance and brings the location information of the obstacle. The numerical implementation details of the iterative scheme is described, and the numerical examples illustrate that the iterative method yields satisfactory reconstructions.
Concerning our future work, the proposed methodology could be extended directly to the case of recovering a sound-hard or impedance obstacle, as well as the three-dimensional case. Although our numerical results show that the proposed novel iterative scheme works very well, the corresponding theoretical justifications of convergence and stability are still open. In other words, the mathematical analysis of this method is beyond the scope of our current work and deserves future investigations. In addition, we would also like to study the applicability of this approach to imaging crack-like scatterers from phaseless data. Moreover, we believe that the reference ball based iteration approach is also a feasible technique for solving the phaseless inverse electromagnetic scattering problem.